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1 .   INTRODUCTION 
1.1.   Scope  of  Thesis 

Multi-level,  or  multi-grid,  methods  have  existed  for  more  than 
fifteen  years,  but  have  only  recently  become  popular.   Thus,  there 
remain  many  unanswered  theoretical  questions.   Several  of  these  are 
addressed  in  this  thesis,  the  major  one  being  the  computational  complexity 
of  multi-level  methods  for  locally  refined  finite  element  grids.   These  are 
grids  for  which  the  ratio  of  diameters  of  the  largest  and  smallest 
elements  is  very  large,  becoming  unbounded  on  the  family  of  discretizations 
of  interest.   The  principal  result  of  this  thesis  is  that  multi-level 
methods  can  approximately  solve  the  linear  algebraic  systems  for  such 
finite  element  discretizations  in  0(N)  operations  for  a  grid  with  N 
elements.   Such  asymptotically  optimal  complexity  bounds  have  been  given 
previously  for  multi-level  methods  on  non-locally-refined  grids,  but  not 
for  locally  refined  grids.   This  result  provides  theoretical  justification 
for  the  use  of  multi-level  methods  on  locally  refined  grids,  which  is 
becoming  increasingly  popular.   The  proof  of  this  result,  which  occupies 
almost  half  of  this  thesis,  is  also  of  interest  for  several  reasons. 
Perhaps  the  most  significant  of  these  is  that  the  rate  of  convergence 
shown  for  the  multi-level  iteration  used  depends  only  on  local  properties 
of  the  finite  element  spaces  involved.   Global  properties,  including 
elliptic  regularity,  are  never  invoked.   As  a  consequence  this  result 
provides  an  independent  verification  of  the  0(N)  complexity  of  multi-level 
methods  for  non-convex  domains,  established  previously  by  Bank  and  Dupont 
(1978).   It  also  opens  the  door  to  easily  computable  bounds  on  the 
convergence  rate  of  multi-level  methods  on  irregular  finite  element  grids. 


In  all,  three  multi-level  algorithms  are  considered  in  detail 
in  this  thesis.   The  first,  analyzed  in  chapter  three,  is  similar  to 
methods  considered  by  Nicolaides  (1977)  and  by  Bank  and  Dupont  (1978) , 
but  is  simpler  than  these  others.   It  is  shown  to  converge  in  0(N) 
operations  for  N-point  quasi-uniform  finite  element  grids  on  convex 
domains,  through  a  convergence  proof  similar  to  those  given  by  the  above 
authors.   This  proof  is  given  here  since  it  is  less  complex  than  an 
earlier  one  based  on  Fourier  analysis  of  the  truncation  error.   The 
algorithm  for  quasi-uniform  grids  is  treated  here  partly  as  background 
for  the  analysis  of  algorithms  for  locally  refined  grids  in  chapter  four. 
It  is  also  of  independent  interest  that  an  0(N)  bound  can  be  shown  for 
this  algorithm  even  though  it  is  simpler  than  those  in  the  literature. 

Two  algorithms  are  considered  in  chapter  four,  both  applicable 
to  locally  refined  grids  on  arbitrary  polygonal  domains.   The  first  of 
these  is  quite  similar  to  the  method  used  in  the  adaptive  finite  element 
code  of  Bank  and  Sherman  (1978).   It  is  the  closest  algorithm  to  theirs 
that  could  be  analyzed  by  the  theoretical  techniques  here.   This  method 
is  shown  to  converge  in  0(N)  operations  subject  to  a  condition  on  the 
rate  of  growth  of  the  dimensions  of  the  finite  element  spaces  used.   The 
necessary  growth  condition  requires  that  the  finer  finite  element  spaces 
used  by  the  multi-level  algorithm  have  far  more  elements  than  the 
coarser  spaces  used.   Such  a  condition  is  also  imposed  in  the  Bank  and 
Sherman  code  though  it  is  less  stringent  there  due  to  their  technique 
of  "level  compression,"  Bank  and  Sherman  (1978). 

Since  these  growth  conditions  can  be  quite  restrictive,  a 
second  algorithm  is  also  considered  in  chapter  four.   It  is  a  modification 


of  the  first,  and  it  will  be  shown  to  converge  in  0(N)  operations  under 
a  considerably  weaker  growth  condition  on  the  dimensions  of  the  finite 
element  spaces  used.   Though  the  theoretical  interest  in  this 
modification  is  apparent,  the  author  is  at. present  unsure  of  the 
practical  utility  of  this  second  algorithm.   One  pays  a  large  penalty 
for  the  weaker  growth  conditions  through  an  increase  in  the  hidden 
constant  in  the  0(N)  work  bound. 

1.2.   History  of  Multi-Level  Methods 

Multi-level  or  multi-grid  methods  are  usually  thought  of  as 
originating  in  the  work  of  Fedorenko  (1961,  1964),  although  related  ideas 
had  been  put  forth  earlier.   The  idea  of  using  several  discretizations 
in  the  solution  of  a  PDE  (partial  differential  equation)  is  actually 
quite  old.   Southwell  (1940)  suggested  using  the  solution  of  an  elliptic 
problem  on  a  coarse  finite  difference  grid  to  obtain  a  starting  value 
for  relaxation  on  a  finer  grid.   He  also  recommended  a  technique  called 
"block  relaxation"  which  is  closely  related  to  multi-level  methods. 
In  this  technique  the  person  doing  the  relaxation  calculation  (in  the 
era  before  computers)  would  find  a  region  of  the  finite  difference  grid 
where  the  residuals  were  all  of  the  same  sign.   A  constant  would  then 
be  added  to  the  value  of  the  trial  solution  throughout  this  region  to 
reduce  the  average  residual  on  this  region  as  much  as  possible.   As  a 
consequence  ,  after  this  correction  the  residuals  would  oscillate  in  sign 
on  this  region.   Southwell  observed  that  as  long  as  there  were  no  large 
areas  where  the  residual  did  not  alternate  in  sign,  relaxation  converged 
rapidly.   Of  course,  though  quite  effective,  such  a  heuristic  approach 
was  not  of  much  help  in  the  early  days  of  computers,  so  simpler  methods 
like  SOR  took  over. 


The  idea  of  making  block  corrections  essentially  died  out  until 

Fedorenko  (1961)  observed  that  rather  than  making  corrections  in  isolated 

blocks  as  suggested  by  Southwell,  one  could  make  a  global  correction  by 

solving  a  so-called  coarse  grid  residual  problem.   For  example,  suppose 

that  in  iteratively  solving  Poisson's  equation  on  a  grid  with  mesh  size 

h  one  has  at  the  n-th  iteration  a  trial  vector  u  £  G,  and  residual 

— n    h 

r  £  G,  where  G,  is  the  vector  space  of  mesh  functions  on  the  grid. 
— n    h        h 

If  A  is  the  discrete  Laplacian  being  used, 

h   _h   .   h 

r  =  f  -  A,  u 
— n   —     h  — n 

where  f^  is  the  data  vector.   Now  let  G~,  be  the  vector  space  of  mesh 

functions  on  a  coextensive  coarser  grid  with  mesh  size  2h  and  suppose  we 

have  an  interpolation  mapping 

X2h:   G2h  *  Gh  • 

One  can  ask  for  a  correction  vector  v   £  G_,  such  that  the  improved 

—      zn 

solution 

h      h   xh   2h 

U  J.1  =  u   +  Tou    v 
— n+1   — n    2h  — 

would  have  smaller  residual.   This  idea  has  been  considered  independently 

by  others,  Wachspress  (1974).   It  is  attractive  since  v   lies  in  the 

lower  dimensional  space  G__  and  so  should  be  relatively  easy  to  compute. 

zh 

9Vi 
Fedorenko' s  fundamental  contribution  was  to  note  that  v   could  be 

chosen  as 

/i  n    -,  \  a     2h   T2h  h 

(1.2.1)  A2hv   =  Ih  r^ 

9Vi 
where  I   :  G,  ->  G_   is  an  interpolation  mapping  and  A?,  is  a  discrete 

Poisson  operator  on  G„,  .   Since  (1.2.1)  is  just  a  discretization  of 

Poisson's  equation  with  data  I,   r  , he  realized  that  (1.2.1)  could  be 

h  — n 


approximately  solved  by  relaxation  and  the  recursive  use  of  corrections  on 
a  still  coarser  grid  G,,  . 

Fedorenko's  theoretical  work  showed  that  for  a  model  problem 
multi-level  iteration  reduced  the  convergence  error  by  an  amount  per 
iteration  independent  of  the  mesh  size.   Thus,  reducing  the  convergence 
error  to  a  magnitude  comparable  to  the  truncation  error  required 
0(N  log  N)  operations  on  an  N  point  grid,  the  same  as  ADI.   Shortly 
thereafter  Bakhvalov  (1966)  considered  beginning  the  solution  process 
on  a  very  coarse  grid,  using  the  solution  there  as  a  starting  value  for 
multi-level  iteration  on  a  slightly  finer  grid,  and  so  on,  gradually 
descending  down  to  the  level  of  refinement  desired.   This  idea, which 
was  probably  old  when  Southwell  was  young,  changes  the  work  bound  from 
0(N  log  N)  to  0(N).   No  other  current  method,  either  direct  or  iterative, 
achieves  this  asymptotically  optimal  work  bound. 

Besides  improving  the  asymptotic  work  bound  on  multi-level 
iteration,  Bakhvalov  extended  the  method  to  a  large  class  of  PDEs  on 
a  square  finite  difference  grid.   After  his  work,  almost  no  work  seems 
to  have  been  done  on  multi-level  methods  until  this  decade.   Then  several 
people  became  interested  in  this  method  including  Brandt,  Hackbusch  and 
Nicolaides.   Brandt  in  particular  is  largely  responsible  for  popularizing 
this  method  and  developing  it  into  a  practical  computational  tool.   The 
other  two  people  mentioned  have  been  working  more  on  the  theoretical 
questions  involved  in  this  method  and  on  generalization  to  finite  element 
grids. 

In  the  last  few  years  there  has  been  a  rapid  growth  of  interest 
in  multi-level  methods,  and  there  is  considerable  ongoing  research.   This 


research  is  directed  both  at  the  theoretical  problems,  Bank  and  Dupont 
(1978),  Hackbusch  (1978),  Nicolaides  (1977),  and  at  developing  algorithms 
for  the  many  applications  areas.   Representative  work  of  the  latter 
type  may  be  found  in  Bank  and  Sherman  (1978),  Brandt  (1977),  Brandt, 
Dendy  and  Ruppel  (1978),  Nicolaides  (1975),  and  Poling  (1977).   If 
the  author  is  not  mistaken,  we  should  see  continued  research  in  this 
area  for  many  years  and  may  eventually  see  multi-level  methods  as 
the  method  of  choice  for  most  parabolic  and  elliptic  problems. 


1.3.   The  Finite  Element  Approach 

Design  and  analysis  of  multi-level  methods  can  be  done  from 
either  a  finite  difference  or  finite  element  point  of  view.   In  one 
sense  these  two  approaches  are  almost  the  same;  finite  element  methods 
are  a  special  class  of  finite  difference  method.   Nevertheless,  the  way 
people  think  about  finite  element  methods  and  the  kinds  of  grids  they 
use  for  these  methods  are  so  different  from  the  usual  way  people 
understand  and  use  finite  difference  methods  that  we  may  think  of  these 
classes  as  completely  separate.   This  thesis  is  devoted  to  analysis 
of  multi-level  methods  from  the  finite  element  point  of  view.   Several 
of  the  many  reasons  the  author  had  for  adopting  this  point  of  view  will 
be  described  in  this  section. 


While  the  multi-level  method  originated  as  a  means  of  solving 
finite  difference  equations  and  was,  until  recently,  solidly  wedded  to 
the  finite  difference  approach,  much  of  the  current  work  on  multi-level 
methods  is  being  done  in  the  finite  element  context.   Starting  with 
Hackbusch  (1977)  and  Nicolaides  (1977) ,  who  seem  to  have  been  the  first 
to  consider  multi-level  methods  as  a  means  of  solving  finite  element 
equations,  there  has  been  increasing  interest  in  this  approach.   The 
reasons  for  this  are  not  complicated.   Finite  element  theory  is  in  many 
respects  superior  to  finite  difference  theory.   A  great  many  problems 
in  the  finite  difference  approach  either  do  not  arise  in  the  finite 
element  approach  or  are  easily  handled  there.   Questions  of  stability, 
the  use  of  irregular  grids,  treatment  of  corner  singularities,  treatment 
of  curved  boundaries  and  construction  of  error  estimates  are  all  much 
simpler  in  the  finite  element  context. 

Aside  from  the  obvious  (to  the  author)  superiority  of  the 
finite  element  method,  there  are  special  reasons  unique  to  multi-level 
algorithms  for  doing  their  theory  in  the  finite  element  setting.   The 
first  of  these  is  that  the  finite  element  error  estimates  are  based  on 
approximation  results  rather  than  on  asymptotic  expansions.   What  this 
means  for  multi-level  methods  is  that,  unlike  some  of  the  early  finite 
difference  arguments  where  rapid  convergence  of  the  iteration  could  be 
shown  only  for  sufficiently  small  size,  no  such  proviso  is  ever 
needed  in  the  finite  element  context.   A  second,  more  important  reason 
is  that  multi-level  algorithms  always  involve  more  than  one 
discretization.   Thus  one  requires  mappings  from  the  space  of  grid 
functions  or  finite  element  functions  of  each  discretization  to  those  of 
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other  discretizations.   In  the  finite  difference  context  such  mappings 
are  done  by  interpolations  which  require  separate  error  analyses.   In 
the  finite  element  context  the  grids  can  be  designed  so  that  the 
finite  element  space  on  a  coarser  grid  will  be  a  subspace  of  the  finite 
element  spaces  on  finer  grids.   Then  the  mapping  from  a  coarser  grid  to 
a  finer  grid  will  be  just  the  natural  injection.   To  go  the  other  way, 
some  kind  of  projection  is  needed,  but  many  natural  ones  are  available. 
The  fact  that  different  finite  element  discretizations  can  be  related 
in  this  way  yields  a  great  simplification  of  the  theory  of  multi-level 
methods.   Unfortunately  only  the  theory  is  simplified,  not  the 
programming.   Rather  than  an  interpolation  formula,  one  has  a  similar 
formula  relating  the  basis  functions  of  one  finite  element  space  to 
those  of  another.   It  is  interesting  that  in  many  cases  these  change 
of  basis  formulas  arising  in  the  finite  element  setting  are  exactly 
the  same  as  the  interpolation  formulas  already  in  use  in  the  finite 
difference  setting. 

There  is  another  reason  why  studying  multi-level  convergence 
is  important  for  this  thesis.   In  chapter  four  we  treat  multi-level 
methods  for  locally  refined  grids.   While  it  may  be  possible  to  analyze 
multi-level  methods  for  locally  refined  finite  difference  grids  in 
special  cases,  it  seems  unlikely  that  an  analysis  as  general  as  that 
given  here  would  be  possible  outside  finite  element  theory.   For 
locally  refined  grids,  finite  element  theory  also  helps  in  coding. 
One  never  needs  to  choose  special  formulas  for  points  where  the  mesh 
size  changes  abruptly  or  worry  about  truncation  error.   There  are 
also  beginning  to  be  reliable  a  posteriori  error  estimates  for  finite 


elements  enabling  one  to  design  robust  adaptive  algorithms,  Babushka 
and  Rheinboldt  (1978).   In  short,  multi-level  methods,  adaptive 
refinement,  and  finite  element  theory  mesh  together  quite  nicely. 

The  main  limitation  of  the  finite  element  approach  is  that  it 
is  far  more  natural  and  powerful  for  self-adjoint  problems  than  for  the 
general  case.   Neither  finite  difference  nor  finite  element  methods  work 
very  well  for  non-self-adjoint  problems,  but  it  is  frequently  easier 
to  see  what  to  do  on  regular  finite  difference  grids  than  on  general 
finite  element  grids.   While  it  would  have  been  possible  to  extend 
many  of  the  results  of  this  thesis  to  the  non-self-adjoint  case  by 
adding  various  restrictions,  that  was  a  more  ambitious  project  than  the 
author  was  prepared  to  attempt.   If  one  takes  the  results  here  as 
indicative  of  the  kind  of  results  that  could  be  proved  in  the  general 
non-self-adjoint  case,  one  should  not  go  far  wrong. 


10 


2.   PRELIMINARIES 
In  this  chapter  we  collect  notations  and  basic  results  from 
finite  element  theory.   Since  this  material  is  completely  standard,  the 
reader  familiar  with  finite  elements  can  probably  skip  this  chapter 
entirely,  referring  back  to  it  later  as  required.   The  reader  less 
familiar  with  the  finite  element  approach  to  partial  differential 
equations  may  find  this  chapter  helpful. 

2.1.   Elliptic  Equations 

Throughout  this  thesis  we  will  be  concerned  with  solving 
self-adjoint  second  order  boundary-value  problems.   The  problem  of  most 
interest  will  be  the  Dirichlet  problem 
(2.1.1a)  Lu  =  f  on  Q 

(2.1.1b)  u  =  0  on  8ft  , 

2 
where  ft  is  a  bounded  open  domain  in  the  plane  IR  with  piecewise  smooth 

boundary  9ft.   L  may  be  any  self-adjoint  second  order  operator,  for 

example  the  Laplacian,  or  more  generally  any  operator  in  divergence  form 

Lu  =  -V  •  (a  V  u)  +  bu 

where  a,  b  are  C   functions  on  ft.   In  order  for  this  operator  to  be 

elliptic  on  ft,  we  require  that  there  be  constants  a.,  a,  b  such  that 

3.  —   a  —   a 

0  <   b  <  b 

hold  throughout  ft.   We  will  also  be  considering  the  Neumann  problem 

(2.1.2a)  Lu  =  f  on  ft 

(2.1.2b)  u  =  0  on  9ft  , 

n 

where  u  is  the  derivative  of  u  in  the  direction  normal  to  the  boundary. 
n 
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In  this  case  ellipticity  requires  that  b  be  positive  on  some  neighborhood 
in  ft.  Bank  and  Dupont  consider  also  the  singular  Neumann  problem,  b  =  0, 
but  we  exclude  it  here. 

The  finite  element  approach  to  boundary-value  problems  is  based 
in  part  on  the  functional  analytic  theory  of  partial  differential 
equations,  in  which  the  boundedness  or  invertibility  of  a  differential 

operator  is  expressed  in  terms  of  Sobolev  norms.   These  norms  are  a 

2 
natural  generalization  of  the  L  norm,  defined  for  non-negative  integers 

I   by 

||u||    =    (       Z         ||Dau||22)1/2    , 
0<|a|<£  L 

where  a  =    (a   ,   a  )    is  a  multi-index  with  a. ,   a     >   0  and 

j  ot  |    =  a     +  a~ 

Da  =        9 


al       a2 

8x        9y 


This  definition  makes  sense  for  all  functions  on  ft  whose  distribution 

2  °o 

derivatives  up  to  order  I   exist  and  lie  in  L  .   The  completion  of  C  (ft) 

in  ||  •  ||  is  denoted  H  (ft)  and  is  a  Hilbert  space  relative  to  the  inner 

product 

(u,v)   =  /  (   Z    (Dau)(D°v))  . 
ft  0<|a|<£ 

If  &   is  a  negative  integer  ||  •  |L  is  defined  by  duality, 

||u||£=   sup  -(^- 


-i      l|v'U 

veH  (ft) 

Here  and  throughout,  (•,  •)  is  the  L   inner  product  while  ||  •  ||  with  no 

2 
subscript  denotes  the  L  norm.   Sobolev  spaces  can  also  be  defined  for 

arbitrary  real  Z   either  through  the  interpolation  theory  of  Hilbert 
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spaces  or  by  Fourier  transforms.   Details  may  be  found  in  Oden  and  Reddy 

(1976). 

Equations  (2.1.1)  and  (2.1.2)  implicitly  require  the  solution 

u  to  be  twice  dif ferentiable,  which  is  overly  restrictive  for  many  purposes 

Consider  the  problem 

(2.1.3a)  Lu  =  f  on  Q 

(2.1.3b)  u  =  0  on  T 

(2.1.3c)  u  =  0  on  I\  , 

n         2 

where  9S7  =  V     U  r  ,  which  includes  both  (2.1.1)  and  (2.1.2)  as  special 

cases.   Let  H  (fi)  be  the  subspace  of  H  (Q)    satisfying  the  (essential) 

Dirichlet  boundary  condition  (2.1.3b)  but  otherwise  unrestricted.   If 

(2.1.3)  is  satisfied  then 

(Lu,v)  =  (f,v),    v  G  Hg(fi)  . 

The  notation  here  means  that  equality  holds  for  all  test  functions  v  in 
H  (fi) .   Defining  the  bilinear  form 

a(u,v)  =  (Lu,v)  , 

we  can  ask  for  the  element  u  in  H  (fi)  such  that 

E 

(2.1.4)  a(u,v)  =  (f,v),    v6Hj(i])  . 

It  is  only  necessary  that  u  lie  in  H   since  Green's  theorem  can  be  used 

E 

to  shift  one  derivative  from  u  to  v.   A  function  u  satisfying  (2.1.4)  is 
called  a  weak  solution  of  (2.1.3),  and  (2.1.4)  is  called  the  weak  form 
of  the  equation.   The  advantage  here  is  that  the  weak  solution  exists 
quite  generally  and  agrees  with  the  classical  solution  satisfying 
equation  (2.1.3)  pointwise  when  the  latter  exists. 
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The  condition  for  existence  and  uniqueness  of  the  solution  of 

(2.1.4)  is  that  the  problem  be  elliptic.   This  means  that  the  energy 
norm  |||  •  |||  defined  by 

|||u|||2  =  a(u,u) 

for  u  E  H  (ft)  is  equivalent  to  the  first  Sobolev  norm  ||  •  ||  .   That  is, 
Jj  1 

ffjM't  <  INII  <  ^H  .   u^  H^ft)  , 

for  fixed  0    ,  G~  >  0  independent  of  u. 

If  the  problem  (2.1.3)  is  elliptic,  a  simple  calculation  shows 
that  for  &  a   non-negative  integer 

l|u|l£+2  >   c  ||Lu||£    ,  u  e  H£+2(ft)    , 

for  some  c  >  0.   The  principal  results  in  elliptic  theory  go  the  other 

3  I 

way.   For  I   >  -  -~  and  any  f  e  H  (ft)  the  weak  solutions  of  (2.1.1)  or 

(2.1.2)  satisfy  the  regularity  estimate 

(2.1.5)  NI£+2  <  CHfll*, 

for  some  c  >  0,  assuming  that  the  domain  is  smooth  and  the  coefficients 

00 

of  L  are  C  .   For  the  technical  definition  of  smooth  domains  see 
Babushka  and  Aziz  (1972)  or  Oden  and  Reddy  (1976). 

The  regularity  estimate  (2.1.5)  implies  existence,  uniqueness, 
and  that  the  solution  has  two  more  derivatives  in  the  mean  square  sense 
than  its  data.   In  particular,  if  f  is  infinitely  dif f erentiable,  u  will 
be  as  well.   These  results  hold  only  for  smooth  domains.   For  domains 
that  are  only  piecewise  smooth  the  corners  introduce  singularities  that 

limit  the  regularity  of  the  problem.   On  a  convex  polygonal  domain  the 

2  I 

solution  will  lie  in  H  but  not  necessarily  in  H  for  £  >  2.   In  this 

case  we  have  the  regularity  estimate 


14 


(2-1.6)  IMU^IIf"* 

3 
for  -  "z-  <  £  _<  0.   For  non-convex  polygonal  domains  (2.1.6)  will  not  hold 

for  I   =  0  because  the  reentrant  corners,  that  is  corners  with  interior 

angles  exceeding  tt,  introduce  severe  singularities.   In  the  worst  case  of 

3    1 
a  slit  domain,  £  must  be  restricted  to  (-  -~,  -  ■=■).   The  solution  has 

less  than  one  and  a  half  derivatives  in  the  mean  square  sense. 


2.2.   Finite  Element  Spaces 

The  finite  element  discretization  of  (2.1.3)  begins  with  the 

selection  of  a  finite  dimensional  subspace  M  of  H  (ft).   The  discrete 

E 

solution  to  (2.1.3)  is  taken  as  that  function  u  E  M  satisfying  the  finite 

element  equation 

(2.2.1)  a(u,<f>)  =  (f,<j>),  <J)  e  M  . 

By  the  Ritz  theorem,  holding  for  self-adjoint  problems,  u  is  the  function 

in  M  nearest  to  the  true  solution  u  in  energy  norm. 

Finite  element  spaces  are  ordinarily  constructed  by  dividing 
the  domain  into  a  finite  family  of  elements,  usually  triangles  or 
rectangles.   Once  the  elements  are  chosen  M  may  be  taken  as  a  finite 
dimensional  space  of  piecewise  polynomial  functions  that  are  polynomial 
on  each  element.   Though  it  is  not  strictly  necessary  for  multi-level 
methods,  throughout  we  will  require  that  M  be  conforming,  that  is  that 
it  lie  in  IL(ft).   This  translates  into  the  requirement  that  the  functions 
in  M  be  continuous  and  satisfy  the  essential  boundary  condition  (2.1.3b) 
exactly.   We  also  require  that  M  be  of  degree  at  least  one  in  order  to 
obtain  convergence.   For  a  finite  element  space  of  piecewise  polynomial 
functions  the  degree  is  the  integer  k  such  that  any  polynomial  of  degree 
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k  can  be  exactly  represented  on  every  element  away  from  the  boundary. 
Finally  It  is  necessary  to  bound  the  highest  degree  polynomial  on  any 
element  and  require  that  M  have  a  local  basis.   This  is  in  order  that  the 
basis  be  uniformly  linearly  independent.   The  question  of  the  linear 
independence  of  the  basis,  or  conditioning  of  the  Gram  matrix  will  be 
taken  up  in  the  next  section. 

Approximation  questions  are  dealt  with  in  finite  element 
theory  in  terms  of  a  sequence  of  finite  element  spaces  of  increasing 
dimension  IM. /._,»  which  eventually  become  dense  in  the  sense  that  for 

any  neighborhood  N  in  H„(fi)  there  exists  J  such  that 

E 

N  n  M.  +   0  for  j  >  J  . 

The  approximation  properties  of  each  of  these  spaces  M.  depend  on  the 

size  of  the  elements  {e}  of  M.  and  the  element  geometry,  as  well  as 

the  degree  of  the  finite  element  space.   Element  geometry  matters  only 

to  the  extent  that  it  must  not  degenerate  as  j  increases.   For  standard 

finite  elements  it  suffices  that  all  elements  In  M.  be  convex  and  have 

3 

interior  angles  bounded  away  from  0  and  IT  uniformly  in  j  . 

The  size  of  an  element  e  is  measured  by  its  diameter,  d  , 

defined  as 

de  =   sup  ||x-y||  2  , 
x,yee      I 

where  ||»||  „  denotes  the  Euclidean  norm.   The  family  of  spaces  {M.}  is 

I  J 

said  to  be  quasi-uniform  if  there  is  a  corresponding  family  of  parameters 

{h.}  and  a  >  1  such  that 

(2.2.2)  -  h.  <  d  <  ah. 

O     3  -     e  -         j 

for  all  elements  e  of  M.  and  all  j.   If  not,  that  is  if  the  constant  O 


16 


must  depend  on  j,  the  family  of  spaces  {M  }  is  said  to  be  locally  refined, 

While  strictly  speaking  these  concepts  apply  only  to  an  infinite  family 

of  finite  element  spaces,  we  will  talk  about  an  individual  finite  element 

space  being  quasi-uniform  or  locally  refined.   This  usage  is  extremely 

convenient  though  it  should  be  born  in  mind  that  this  really  makes  sense 

only  with  respect  to  membership  in  an  infinite  family  of  spaces {M. }. 

This  definition  of  the  mesh  parameter  h.  associated  with  a 

quasi-uniform  finite  element  space  M.  is  not  completely  standard.   Often 

h.  is  taken  as  the  maximum  of  the  element  diameters.   However,  the 
J 

extra  freedom  allowed  by  (2.2.2)  makes  no  difference  in  the  error 

estimates  since  a  is  a  fixed  constant,  and  is  quite  convenient.   In 

particular  this  freedom  makes  it  easy  to  have  a  quasi-uniform  family 

■f^ili  /n  1N  parameterized  by  h  rather  than  j.   In  this  case  the  diameters 
tl  n£  (,U,  L) 

of  elements  in  M,  are  required  to  satisfy 

—  h  <  d   <  Gh  . 
a   —  e  — 

The  standard  finite  element  error  estimates  give  the  rate 
of  convergence  in  any  Sobolev  norm  of  the  finite  element  approximations 
in  a  quasi-uniform  space  M  in  terms  of  the  mesh  parameter  h  and  the 
degree  of  M  .   For  smooth  solutions  the  convergence  rate  may  be  made 
arbitrarily  high  by  taking  finite  element  spaces  {M,  }  of  sufficiently 
high  degree.   But  on  a  polygonal  domain  the  corner  singularities  limit 
the  rate  of  convergence  attainable  by  any  family  of  piecewise 
polynomial  finite  element  spaces. 

The  most  important  error  estimate  in  the  next  chapter  will  be 
the  L  error  estimate 


17 


(2.2.3)  ||u-uh||2    <   Ch2||f||2, 

L  L 

which  holds  for  quasi-uniform  spaces  {M,  }  of  degree  at  least  one  when 

h 

2 

the  partial  differential  equation  is  H  regular,  that  is,  when  (2.1.6) 

holds  with  I  =   0.   Thus,  this  bound  applies  to  the  Dirichlet  problem 
(2.1.1)  and  the  Neumann  problem  (2.1.2)  on  any  convex  polygonal  domain. 

The  convergence  in  energy  for  these  two  problems  is  of  first 
order  on  any  convex  polygonal  domain.   On  non-convex  domains  it  is  of 
less  than  first  order.   For  a  domain  with  a  slit,  the  worst  case 
ordinarily  arising,  it  is  of  order  one-half. 

Much  of  the  utility  of  the  finite  element  method  comes  from 
the  fact  that  these  meager  convergence  rates  attained  on  polygonal 
domains  are  easily  improved  through  the  addition  of  singular  functions 
or  by  doing  local  refinement  to  cope  with  the  corner  singularities.   In 
chapter  four  the  use  of  multi-level  methods  on  locally  refined  finite 
element  grids  will  be  considered.   A  simple  error  estimate  showing 
the  quality  of  approximation  of  finite  element  methods  on  such  grids 
will  be  given  there. 

2.3.   Linear  Equations 

By  selecting  a  basis  for  the  finite  element  space  M,  the 

finite  element  equation  (2.2.1)  is  translated  into  a  linear  system. 

Let  M  be  a  finite  element  space  for  problem  (2.1.3)  and  let  l<p  /._- 

be  a  basis  for  M.   Any  function  u  6  M  can  be  expressed  in  this  basis  as 

N 
u  =  £   u.  (f>.  , 
i=l  -1   X 
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where  u^  is  the  coordinate  vector  of  u  relative  to  this  basis.   The 

convention  of  using  underbars  for  vectors  will  be  adhered  to  throughout. 

Given  the  basis  {<}>.},  we  can  define  the  mass  and  stiffness 

1 

matrices  M  and  K  with  elements 

mij  =  C*±.  Oj) 

kij  =  a(V  V 

respectively.   M  is  just  the  Gram  matrix  of  the  basis  relative  to  the 

2 
L  inner  product  while  K  is  the  Gram  matrix  in  the  energy  inner 

product.   To  form  the  linear  system  corresponding  to  the  Ritz  equation 

(2.2.1),  let  f  e  L  (fi)  and  let  F  be  given  by 

F±  =    (f,  (JO,     i  =  1,...,  N  . 

Then  (2.2.1)  is  equivalent  to  the  linear  system 

(2.3.1)  Ku  =  F 

where  jl  is  the  coordinate  vector  of  the  finite  element  solution. 

Numerical  solution  of  the  linear  system  (2.3.1)  is  aided  by 
several  properties  the  matrices  K  and  M  have.   First,  M  will  always  be 
symmetric  positive  definite  and  since  L  is  self-adjoint  K  will  be 
symmetric  positive  definite  as  well.   Second,  for  the  usual  finite 
element  bases  M  and  K  will  be  quite  sparse.   Finally,  both  M  and  K  are 
usually  well  conditioned.   The  common  finite  element  bases  are  largely 
motivated  by  these  last  two  considerations  which  are  of  great  importance 
for  the  numerical  inversion  of  (2.3.1)  or  (2.3.2)  whether  this  is  done 
directly  or  iteratively. 

To  make  these  considerations  precise,  consider  a  family  of 

finite  element  spaces  {M . }   -.of  increasing  dimension.   We  require  that 

J  C)  N- 

each  space  M.,  1  <  i  <  °°,  have  a  basis  {(J).   }.  ,  that  is  local  in  the 
j    —  i   i=l 
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sense  that  the  number  of  non-zero  elements  in  each  row  of  the  corresponding 
mass  and  stiffness  matrices  is  bounded  uniformly  in  j .   This  will  be  true 
if  the  number  of  basis  functions  in  {<£>.   }._,  whose  supports  intersect 
the  supports  of  any  given  basis  element  <$.         is  bounded  uniformly  in  j  . 
Usually,  under  these  conditions  the  condition  number  k(M.)  of 
the  mass  matrix  M.  ,  1  j<  j  <  °°,  will  be  uniformly  bounded.   For  this  to 
hold  the  proper  scaling  must  be  chosen  and  the  condition  numbers  of 
the  element  mass  matrices  must  be  uniformly  bounded.   This  is  true  for 
the  common  finite  element  bases  as  long  as  the  element  geometry  does 
not  degenerate.   Then  if  the  basis  elements  are  scaled  as 

||^j)||  =  1,    1  <  i  <  N.  ,    1  <  j  <  °°  , 

which  we  assume  throughout,  it  is  immediate  that  the  eigenvalues  of 
M.  will  lie  in  an  interval  [— ,  o]    for  O  >   1  independent  of  j. 

Since  the  stiffness  matrices  {K. }  are  exactly  analogous  to 
the  matrices  arising  in  finite  difference  discretizations,  their 
condition  numbers  cannot  be  uniformly  bounded.   Just  as  in  the  finite 
difference  case  the  maximum  eigenvalue  of  K.  satisfies 

(2.3.2)  A    (K.)  <  ch"2 

max  j   —   mm 

under  ordinary  circumstances,  where  h  .   is  the  diameter  of  the 

mm 

smallest  element  of  M. .   The  minimum  eigenvalue  of  K.  can  be  bounded 

3  3 

below  by  a  simple  argument.   We  have 

T  T         T 

x  K.x_       x  K.x.     x  K.x 

A  (K.)  =  min  — ^—  >   min    J   min  — ^~ 

J         xx        x  M.x      x  x 
_  -,_      

=  Am:1  K.)  XAK.) 
1  3        J    1  3 

>   A. (MT1  K.)/a  . 
~  1   J    3 
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The  operator  M   K   is  the  linear  operator  corresponding  to  the  finite 
element  equation  (2.2.1).   Since  the  finite  element  eigenvalues 
approximate  the  eigenvalues  of  the  PDE  from  above,  A  (K.)  must  be 
bounded  below  uniformly  in  j .   Thus  one  ordinarily  has 

(2.3.3)  k(K.)  <  ch"2 

2     -       mm 

or   for  a  quasi-uniform  spaces  JVL  . 

(2.3.4)  k(K.)   <   ch"2    . 

The  difficulty  in  iteratively  inverting  the  linear  system  (2.3.1)  arises 
directly  from  this  growth  of  the  condition  of  the  stiffness  matrix  as 
h  decreases. 

These  bounds  on  the  condition  of  M  and  K  are  proved  in 
Fried  (1971)  and  may  be  found  also  in  Strang  and  Fix  (1973). 
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3.   QUASI-UNIFORM  GRIDS 
In  this  chapter  a  multi-level  algorithm  simpler  than  those 
considered  in  previous  theoretical  work  will  be  given.   It  will  be 
shown  to  be  of  optimal  order, producing  a  good  solution  on  a 
quasi-uniform  finite  element  grid  in  0(N)  operations  on  an  N 
point  grid. 

3.1.   Introduction 

In  this  chapter  we  look  at  a  multi-level  algorithm  for 
self-adjoint  elliptic  boundary  value  problems  on  quasi-uniform  finite 
element  grids.   This  algorithm  is  related  to  algorithms  considered  by 
Nicolaides,  Bank  and  Dupont ,  Hackbusch  and  others,  but  is  less  complex 
than  these  other  algorithms.   Like  many  of  those  considered  by  others, 
the  algorithm  here  is  of  optimal  order,  producing  a  second  order 
accurate  solution  in  0(N)  operations  on  an  N  point  grid. 

The  optimal  order  convergence  of  the  algorithm  in  this  chapter 
was  originally  shown  for  finite  difference  discretizations  by  a 
cumbersome  Fourier  technique,  and  then  for  finite  element  discretizations 
in  the  same  way.   Here  an  easier  proof,  following  the  ideas  of  the  above 
authors,  will  be  given.   The  proof  is  based  particularly  on  the  work  of 
Bank  and  Dupont  but  employs  several  minor  simplifications.   Throughout 
we  follow,  for  the  most  part,  their  notation.   It  is  convenient  here 
since  little  reference  to  algebraic  quantities  such  as  vectors  or 
matrices  will  be  made. 
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2 
3.2.   L  Convergence 

In  this  section  and  the  next  we  examine  the  convergence  of  a 

simple  multi-level  algorithm.   This  algorithm  can  be  used  to  obtain 

finite  element  approximations  to  the  solutions  of  the  elliptic 

problems  (2.1.1)  and  (2.1.2)  on  a  convex  polygonal  domain  ft.   The 

restriction  to  convex  domains  is  necessary  since  we  need  to  use  the 

2  2 

L  error  estimate  (2.2.3)  which  requires  H  regularity.   Using  curved 

elements  as  in  Bank  and  Dupont  [1978]  the  convergence  results  given 

here  extend  immediately  to  piecewise  smooth  domains  free  of  reentrant 

corners.   For  details  on  the  treatment  of  curved  elements,  the 

interested  reader  is  referred  to  their  paper. 

The  algorithm  considered  here  uses  a  nested  family  of 

quasi-uniform  finite  element  spaces  iM.}.-,  nested  in  the  sense  that 

M.  C  M.  ,.  ,    1  <  j  <  °o  . 
3  3+1 

To  approximately  solve  the  finite  element  equations  for  one  of  these 
spaces  M  ,  k  _>  2,  this  algorithm  uses  a  simple  simultaneous  displacement 
smoothing  iteration  applied  to  the  equations  for  M,  plus  approximate 
solution  of  related  problems  in  M,  - .   If  k  -  1  >  2  these  approximate 
solutions  are  obtained  recursively  by  applying  the  same  algorithm  to 
the  equations  for  M,  ,.   The  recursion  continues,  finally  requiring 
approximate  solutions  to  the  equations  for  M  .  These  can  be  found 
directly  or  by  some  convenient  iteration. 

The  multi-level  method  considered  here  is  related  to  algorithms 
for  finite  elements  by  Nicolaides,  Hackbusch,  and  Bank  and  Dupont  and  to 
methods  for  finite  difference  equations  considered  by  Brandt,  Hackbusch, 
Nicolaides,  Bakhvalov,  Fedorenko  and  others.   The  algorithm  here  differs 
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from  algorithms  described  by  others  in  that  the  solution  to  the  equations 

for  M,  requires  only  two  approximate  solutions  to  the  equations  for  M,  _,  • 

In  other  multi-level  algorithms,  to  solve  the  equations  for  M, 

approximate  solution  of  related  problems  in  M,  ,    is  used  in  two  ways. 

First  it  is  used  as  a  means  of  obtaining  a  good  starting  value  for 

iteration  on  M  .   Subsequently  it  is  used  as  an  acceleration  procedure 

for  iteration  on  M,  . 
k 

The  first  approximate  solution  of  the  equations  for  M,  - 

is  used  to  obtain  a  starting  value  for  iteration  on  M,  ,  and  considerable 

accuracy  is  required.   It  is  necessary  to  reduce  the  convergence  error 

of  the  problem  on  M    so  it  is  comparable  to  the  truncation  error  on 

M,  n .   In  subsequent  solutions  of  related  problems  in  M,  ,  used  to 
k-1  r  k-1 

accelerate  convergence  of  the  iteration  in  M  less  accuracy  is  required. 

K. 

For  example,  in  one  of  the  proofs  given  by  Bank  and  Dupont  and  in  the 
related  proof  of  theorem  4.1  in  chapter  four  of  this  thesis  it  suffices 
to  solve  the  related  problem  with  a  relative  accuracy  of  33%. 

The  point  of  this  chapter  is  that  if  all  of  the  related 
problems  in  M    are  solved  to  the  same  accuracy,  namely  the  error  in 
the  approximate  solutions  must  be  comparable  to  the  truncation  error 
in  M,  . ,  then  only  two  approximate  solutions  in  M    are  required. 
To  approximately  solve  the  elliptic  problem  (2.1.1)  or  (2.1.2)  in 
M,  ,  k  _>  2,  the  algorithm  begins  by  obtaining  a  good  approximate 
solution  to  the  finite  element  equations  for  M   .   This  approximate 
solution  is  used  as  a  starting  value  for  a  smoothing  iteration  on  M  . 

K. 

Though  this  iteration  could  be  continued  until  adequate  convergence 
was  attained  this  would  be  inefficient.   Instead,  after  a  fixed  number 
of  iterations,  m,  independent  of  k,  we  terminate  the  iteration. 
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The  resulting  approximate  solution  is  then  corrected  by  subtracting  from 
it  the  projection  of  the  convergence  error  onto  M   .   This  corrected 
approximate  solution  in  M  is  then  taken  as  the  accepted  solution. 

The  reason  this  is  effective  goes  back  to  the  work  of 
Southwell.   He  observed  that  at  least  initially  relaxation  rapidly 
reduces  the  residual,  while  the  error  itself  may  be  only  slowly  reduced. 
At  the  end  of  our  algorithm's  smoothing  iteration  the  residual  will  be 
small  relative  to  the  convergence  error.   This  means  that  the 
convergence  error  is  quite  smooth.   In  this  case  the  convergence  error 
can  be  well  approximated  in  the  coarser  space  M    and  the  final 
correction  will  remove  most  of  it. 

Brandt  describes  this  same  phenomenon  in  terms  of  the  reduction 
of  Fourier  components  of  the  convergence  error.   Most  convergence 
proofs,  including  ours,  follow  Fedorenko  who  considered  the  reduction 
of  the  error  components  in  the  direction  of  eigenfunctions  of  the 
discrete  partial  differential  operator.   Fedorenko  refers  to  "good" 
meaning  smooth  eigenfunctions  and  "bad"  meaning  oscillatory  eigenfunctions, 

For  practical  purposes  many  possible  smoothing  iterations 
would  be  acceptable.   Jacobi,  Gauss-Seidel,  and  symmetric  SOR  are  most 
often  used.   The  author  expects  that  improvement  could  be  made  through 
the  use  of  Chebyshev  acceleration,  although  in  numerical  tests  so  far 
changing  to  more  subtle  iterative  schemes  has  had  little  effect  for 
the  well  behaved  self-adjoint  problems  considered  here  (Bank  and 
Sherman  [1979],  Brandt  [1977]). 

In  this  section  the  smoothing  iteration  will  be  taken  as  a 
simple  but  computationally  unattractive  simultaneous  displacement 
iteration  considered  also  by  Bank  and  Dupont.   The  advantage  of  this 
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iteration  is  that  the  resulting  convergence  is  easy  to  analyze  and  is 

basis  independent.   In  the  next  section  we  will  show  that  this  smoothing 

iteration  can  be  replaced  by  a  more  computationally  attractive 

iteration  without  altering  the  convergence  result.   In  this  case  the 

algorithm  will  be  shown  to  be  of  optimal  order,  producing  a  second 

order  accurate  solution  in  0(N)  operations  on  a  finite  element  grid 

with  N  elements. 

Let  (M . } .  ,  be  a  nested  family  of  quasi-uniform  finite 
J  3  =  1 

element  spaces,  and  suppose  we  have  associated  mesh  parameters  ih,} 
satisfying 

(3.2.1)  h.  =  p1"3!^  ,     j  >  1  , 

for  p  >  1  a  fixed  constant.   Let  u    G  M,,  j  >  1,  be  corresponding 
finite  element  solutions  of  problem  (2.1.1)  or  (2.1.2).   That  is, 
u    G  M.  are  functions  satisfying 

(3.2.2)  a(u(j),(b)  =  (f,<J))  ,    <f>  e  M   . 
From  the  error  estimate  (2.2.3) 

(3.2.3)  ||u-u(j)||<  chJllfH 

for  fixed  c  >  0  where  u  is  the  true  solution.   In  the  space  M.  associated 

3 

with  the  finite  element  equation  (3.2.2)  we  have  eigenf unctions 

(')  N* 
{i> .      }  ^  ,  where  N.  is  the  dimension  of  M.,  and  associated  eigenvalues 

i    i^l         3  J 

U^}.3   with 
i   i=l 

0  <  A<j)  <  A<j)  <...<  X(Tj)  . 
1   —   2   —   —  N. 

3 

That  is 
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a(^j),<}>)  ■  A<j)(^j),4>)  ,    <fr  e  M  , 


for  1  £  i  <_  N..   From  the  eigenvalue  bound  (2.3.2)  and  the  fact  that 
the  eigenvalues  of  M  ,  the  mass  matrix  corresponding  to  M  ,  lie  in 
[— ,  a]  for  a  >  1  independent  of  j,  we  have  the  bound 

(3.2.4)  AX7  <  c  hT2  . 

To  solve  the  finite  element  equation  in  M,  ,  k  >^  2,  the 

algorithm  to  be  considered  produces  a  sequence  of  iterates  u  ,  u9,..., 

u    in  M  beginning  with  an  approximate  solution  u_  to  the  finite 
nrrJL     k.  U 

element  equations  in  M,  - .   For  some  fixed  positive  integer  m,  the 

approximate  solution  u  in  will  be  taken  as  the  final  answer.   The 

m+1 

iterates  u.  ,  1  <  I   _<  m,  are  produced  by  a  simple  smoothing  iteration 
while  the  last  one,  which  will  be  the  accepted  solution,  is  the  result 

of  a  correction  involving  the  approximate  solution  of  a  related 

2 
problem  in  M,  , .   For  a  fixed  constant  c_  >  0  and  data  f  in  L  (Q)    the 
k-1  0 

algorithm  is  as  follows: 


1.   Let  u  G  M    be  an  approximate  solution  to  (3.2.2) 


satisfying 


(k-1)  I,  . 
u  -u       < 


(3.2.5)  HVU   Hi  c0  (hk-l}  llfll 

(k-1) 
where  u      is  the  solution  of  (3.2.2). 

2.  Let  u„ ,  1  <_  £  <_  m  be  the  element  of  M  satisfying 

(3.2.6)  (u£-u£_lf<J))  =  -^y  (a(u£_1,(t.)  -  (f,(J)))  ,     *  e  ^  • 

X 

3.  Let  v  £  M,  ,  be  the  solution  of 

k-1 

a(v,<£)  =  a(u  ,<J))  -  (f  ,<J>)  ,    <J>  G  M    . 
m  k-1 
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That    is ,    v   is    given  by 

a(v,cf>)    =   -(r,(f>)    ,  <j>  £  Mk_1 


where,  r  £  M,  is  defined  by 
k 

(r,<|>)  =  (f ,<)))  -  a(u  ,({)),    f£M,  . 

in  k. 

Let  v  be  an  approximation  to  v  satisfying 

(3.2.7)  ||v-v||  <   cQ  h^  ||r|| 

and  let  u  , ..  be  given  by 
m+I 

(3.2.8)  u  ..  =  u  -  v 

m+1    m 

The  corrected  approximate  solution  u  , .  is  taken  as  the  final  answer. 

m+1 

If  we  can  show  that  the  approximate  solution,  u   . ,  produced 
is  accurate  enough  uniformly  in  k  >^  2,  the  use  of  this  algorithm 
recursively  to  solve  the  subproblems  in  steps  one  and  three  will  be 
immediately  justified.   This  is  the  content  of  the  following  theorem. 

Theorem  3.1.   For  any  cn  >  0  there  is  a  positive  integer  m, 

independent  of  k,  such  that  the  approximate  solution  u  .,  produced  by 

m+1 

algorithm  (3. 2. 5)-(3. 2. 8)  satisfies 

(3.2.9)  H«1|H-.<k>ll<c0bJ||f||. 

Proof.   For  convenience  we  drop  the  superscript  k  on 
eigenfunctions  and  eigenvalues  and  the  dimension  N,  of  M,  .   From  (3.2.3) 
and  the  triangle  inequality 

||u(k)   -.^"llicCtaJ  +  bj^Hfll 

<    c(l   +   p2)    h2.  ||f||    . 
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From  this   and    (3.2.5) 

(3.2.10)  ||u0  -  u(k)||  <    (c(l  +  p2)   +  cQ   p2)    h2  ||f||   . 

Now,    for  ft  =   0,    1,...,    m+1   let   ep   be   the   convergence   error 

(k) 

e£  =  Ui  ~   U  * 

Expanding  e_  in  eigenfunctions 

N 


0   1=1  i  i 


e^  =  E  a..  \p_. 
Then  by  (3.2.6) 


N  X. 

e  =  Z  a.  (1  -  t— )  \i>. 

i=l         N 

It  remains  to  estimate  the  result  of  the  correction  in  step  3, 

We  have 

N  A. 

r  =  E  a.  A.  (1  -  -r— )  \b , 

.-IX      A     Tl 

1=1  N 

and  by  the  error  estimate    (2.2.3) 

llv-eJIich^Hrll; 


so  by    (3.2.7) 


Vill-HVlli  »•.-'« +  »«i 


<ch2_1||r||+c0h^1||ri 
<    (c  +  cQ)   p2  h2.  ||r|| 


=   (c  +  c0)   P2<|"     a.   A.C1-    %**. 

1=1  N 


<    (c  +  O    p2  h2  ||en||       max      (A(l  -  -^)m) 
U  k       U     0<A<AXT  AN 
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2 
where  the  last  step  depends  on  the  L  orthogonality  of  the  eigenf unctions 


{\p  .}  .   By  a  simple  calculation 


1 


max      (A(l  -  -— )    )  _<  — -   . 

0<X<XXT  AN  "  m¥1 

N 


Then  from  the  eigenvalue  bound  (3.2.4) 


II     II*    2  C  +  CQ  ||   „ 

II Villi  c  p   i5a~lleoll  ' 

It  follows  from  this  and  (3.2.10)  that 

llVilliep2     ^      (c(l  +  p2)+c0p2)||f||. 

Choose  m  large  enough  that 

2 

(3.2.11)       m+1  >  -2-2-  (c  +  c_)  (c(l  +  p2)  +  c.  p2)  , 

-  cQ        0  0 

and  the  result  follows.  □ 

In  this  proof  we  have  carried  the  constants  p  and  c  explicitly 

to  show  the  dependence  of  m  on  these  quantities.   From  (3.2.11)  m 

4 
apparently  goes  as  —  .   This  sharp  dependence  of  the  number  of  iterations 

C0 
m  on  the  mesh  ratio  is  typical  of  multi-level  algorithms,  but  the 

dependence  on  c   is  not.   For  more  common  multi-level  algorithms,  such 

as  that  considered  by  Nicolaides  [1977]  and  that  considered  in  chapter  4 

for  locally  refined  grids,  m  depends  on  c   only  as  |log(cn)|.   Thus  if 

one  wishes  to  reduce  the  convergence  error  well  below  the  truncation 

error,  the  algorithm  here  cannot  be  recommended.   However,  in  the  usual 

case  where  it  suffices  to  have  convergence  error  comparable  to  the 

truncation  error  this  algorithm  should  perform  about  as  well  as  other 

multi-level  methods. 
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3.3.   Computational  Cost 

While  iteration  (3.2.6)  used  in  step  two  of  algorithm  (3.2.5)- 

(3.2.8)  was  convenient  in  the  last  section  since  it  made  the  analysis 

there  basis  independent,  it  is  impractical.   Each  of  the  iterations, 

(3.2.6),  required  the  solution  of  a  linear  system  involving  the  mass 

matrix  because  of  the  implicit  definition  on  the  left  side.   To  avoid 

this,  we  will  show  in  this  section  that  the  smoothing  iteration  (3.2.6) 

can  be  replaced  by  a  simple  simultaneous  displacement  iteration  without 

impairing  the  convergence  result.   This  analysis  follows  Bank  and 

Dupont  [1978].   With  this  replacement  the  overall  algorithm  is  of  optimal 

order,  producing  a  second  order  accurate  solution  in  0(N)  operations  for 

a  finite  element  grid  with  N  elements. 

Given  a  basis  {(J).}.  ,  for  M,  ,  for  each  u  £  M,  there  corresponds 
l  i=l      k  k 

N 
a  coordinate  vector  u  G  |R  .   As  in  the  last  section,  we  have  temporarily 

dropped  the  subscript  k  on  the  dimension  N  of  M  and  on  eigenvectors  and 

eigenvalues.   Let  b(»,  •)  denote  the  bilinear  form 

T 
b  (u,v)  =  u.  v 

The  norm  induced  by  this  bilinear  form  is  the  same  as  the  Euclidean  norm 

on  corresponding  vectors, 

11-11?   -  b(u,u)    =  ||u||2      . 

I 

Now  let  {_$.}  be  the  eigenvectors  of  the  stiffness  matrix  K  with 


corresponding  eigenvalues 


0<  X  <  x  «...«  xN 


The  eigenvectors  {^.}  are  coordinate  vectors  for  eigenfunctions  {i|>.} 
satisfying 
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(3.3.1)  a(ifi.,    $)    =   A.    b(^.,    <j>)    ,  <J>  G  M 

XXX  K. 

and   the  orthogonality  relations 

a($   ,   $.)   =  bC^,   J.)  -  0  for  i  ^  j    . 

We  wish  to  consider  algorithm  (3. 2.5)-(3. 2, 8)  of  the  last 
section  but  with  the  smoothing  iteration  (3.2.6)  replaced  by  the 
simultaneous  displacement  iteration 

(3.3.2)  b(u£  -  u^_1,  ((.)  =  -±-   (a(u£_1,  cf>)  -  (f,  $))  ,    <j)  S  MR  . 

The  use  of  X„  here  is  only  for  convenience.   It  could  be 

N 

replaced  by  a  convenient  upperbound  on  the  eigenvalues  of  the  stiffness 
matrix,  such  as  (2.3.2),  without  altering  the  convergence  result. 
The  convergence  of  this  algorithm  is  shown  by 

Theorem  3.2.   For  any  c   >  0  there  is  a  positive  integer 

(k) 
m  independent  of  k  >  2  such  that  the  approximate  solution  u  * 

—  m+1 

produced  by  algorithms    (3. 2. 5  (- (3. 2. 8)   with    (3.2.6)    replaced  by    (3.2.2) 
satisfies 

llVl-u«||<cohj2||£||. 

Proof.      The   proof    is   a   slight   modification   of   the   proof   of 
Theorem   3.1.      Expanding   e      in   the   eigenfunctions    {^.}    instead   of    in   the 
eigenfunctions   {^.}  we  have 

N 

e_   =      E      a.    $ .    . 
0        .    ,      x     l 
i=l 


Then,    as   before, 


N  \. 

e     =      E      a .  (1  -  - — )      y .    . 

m       i=l     x  An 

k 
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Let  r  G  M  be  given  by 


and  r  G  M.  by 
k 


Then 


so  we  have 


(r,<J>)  =  -  a(e  ,<}>),    f  £  M, 
m  k 


b(r,<}>)  =  -  a(em,  <j>)  ,    ♦  G  M   • 


r  =  Mr  ; 


2     T         — 1   T      1 
rll  =  r  Mr  =  (M  r)   M(M  f) 


~T  -1. 
=  r  M  r 


««llllP2-«ll«l6 


As  in  the  proof  of  Theorem  3.1 


Now 


Vi'l-  (c(1  +  p2)  +  co  p2)  hkl|r|1 

<  (C(l  +  p2)  +  CQ  p2)a1/2  h2  Hr^ 


1=1  AN 

<||eJL   max    (X(l  -  ^-)m) 

0<A<X„        AN 

N 


N 


where  the  first  inequality  uses  orthogonality  of  the  eigenfunctions  {^.} 
in  the  inner  product  b(*,  •)•   Then  since 
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,1       .  .2  T  1/2      .T        1/2      v 

lle0ll     =   *q  %0   =    (  -^0)       (M       ■^0) 

(3.3.3) 

IK,1/2         i|2  1  ||       ||2  1  n       |,2 

=  I|M  ^11  2    >  -  N-qII  2   =  ^  He0l!b 

A/  A/ 

and  since  the  maximum  eigenvalue  L,  satisfies  the  bound  (2.3.2)  the  rest 

N 

of  the  proof  goes  through  exactly  as  before.  □ 

This  theorem  shows  that  by  twice  approximately  solving  the  lower 

dimensional  linear  system  for  M  ..  and  doing  a  fixed  number  of  simultaneous 

displacement  iterations  on  M.  one  can  obtain  a  second  order  accurate 

k 

solution  in  M  .   Since  the  convergence  shown  by  Theorem  3.2  is  independent 

K. 

of  k  >_  2,  the  two  related  problems  in  M    can  be  solved  by  recursively 
applying  the  same  algorithm  to  their  smaller  linear  systems.   Continuing 
the  recursion  one  ends  up  repeatedly  having  to  solve  the  linear  system  for 
JVL  .   This  can  be  done  either  directly  or  by  any  convenient  iterative  method. 

To  estimate  the  computational  cost  of  this  algorithm  observe 
that  from  (3.2.1)  the  dimensions  of  the  finite  element  spaces  {M.}  must 
satisfy 

C,  p  **  1  N  1  c2  p  J  ,     1  £  j  <  °°  , 

for  c  ,  c_  >  0  independent  of  j.   This  assumes  that  h  ,  the  mesh  size  of 
M  ,  is  fixed  and  the  degree  of  the  highest  degree  polynomial  on  any  element 
is  bounded  uniformly  in  j.   Then,  if  s_  is  the  storage  required  for  the 
direct  or  iterative  inversion  of  the  linear  system  for  M  ,  the  total 
storage  must  be  proportional  to 

k  k 

s  +  Z   N.  <  sn  +  c    E   p  3 
0    .  0  3   —     0    2  .  - 
3=2  J  J=2 

2k 

<  s„  +  c, 


-  °     2  i-p-2 
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for  p  >  1.   Since  s   is  fixed,  assuming  the  mesh  ratio  p  is  greater  than 
one,  the  total  storage  required  depends  only  linearly  on  the  dimension  of 

the  solution  space  M.  . 

k 

To  bound  the  computation  time,  observe  that  the  time  required 
to  do  one  of  the  simultaneous  displacement  iterations  (3.3.2)  for  M.  is 
proportional  to  N..  Let  t  be  the  computation  time  for  the  approximate 
solution  of  the  linear  system  for  M  .  Then  taking  into  account  that  we 
must  recursively  apply  the  algorithm  twice  to  M,  - ,  four  times  to  M ,_9» 
and  so  on  the  computation  time  T  is  bounded  by 

V  k    lr  ' 

(3.3.4)  T  <  2  tA  +  cm  I     2  J  N. 

"     °      3=2       3 

<  2k  t  +  cc  m  E   2k"j  p2k 
j-2 

^k  2k  „   ,  2  vk--j 

£  2   t  +  cc2  mp    E   (-j)   J  . 

j=2  P 

We  have  neglected  the  computational  cost  of  carrying  out  the 
mappings  of  functions  on  one  grid  to  another  in  steps  one  and  three  of 
algorithm  (3. 2. 5)- (3. 2. 8) .   Despite  the  fact  that  the  spaces  {M.}  are 
nested,  the  mappings  required  in  steps  one  and  three  are  not  free  since 
in  general  each  space  M.  uses  a  different  basis.   It  is  however  easy  to 
see  that  the  above  bound  on  the  computation  time  is  unaltered  by  taking 
into  account  this  additional  cost. 

Now  from  (3.3.4)  we  have  the  bounds 

0(N,  )  for  p  >  /2 

k 


T  < 


0(N  log(N  ))    for  p  =  /2 


log  2 


.  0  N/  XOg  p      for  1  <  p  <  /2  . 
^   k 
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The  natural  choice  p  =  2  gives  an  optimal  order  algorithm  and  is  convenient 
in  programming. 
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4.   LOCALLY  REFINED  GRIDS 
This  chapter,  which  contains  the  principal  results  of  this 
thesis,  presents  two  multi-level  algorithms  achieving  0(N)  complexity 
bounds  on  locally  refined  grids.   One  of  these  attains  this  optimal 
complexity  bound  under  weaker  than  expected  restrictions  on  the 
dimensions  of  the  finite  element  spaces  used  by  the  method.   These 
optimal  order  convergence  results  are  a  consequence  of  an  approximation 
theorem  proved  in  this  chapter.   Aside  from  its  use  in  extending  multi- 
level convergence  theory  to  locally  refined  grids,  this  approximation 
result  is  of  interest  since  it  relies  only  on  local  properties  of  the 
finite  element  spaces  used  and  is  not  based  on  elliptic  regularity. 
As  a  consequence  it  provides  an  independent  verification  of  the  0(N) 
convergence  of  multi-level  methods  on  non-convex  domains  shown  previously 
by  Bank  and  Dupont . 

4.1.   Introduction 

There  are  a  variety  of  reasons  why  one  might  wish  to  use  locally 
refined  grids.   They  can  be  used  to  maintain  high  order  accuracy  in  the 
face  of  singularities  caused  by  discontinuous  coefficients  or  the  corners 
of  the  domain,  or  to  resolve  boundary  layers.   They  may  simply  result 
from  an  adaptive  discretization.   They  can  also  be  used  to  give  good 
resolution  of  the  part  of  the  solution  of  greatest  interest  without 
incurring  the  computational  cost  of  a  fine  global  grid.   This  situation 
occurs,  for  example,  on  unbounded  domains.   As  long  as  the  ratio  of  the 
diameter  of  the  largest  element  to  that  of  the  smallest  remains  uniformly 
bounded  on  the  family  of  finite  element  spaces  of  interest,  multi-level 
convergence  results  for  quasi-uniform  grids,  such  as  those  in  the  last 
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chapter,  remain  valid.   However,  in  many  cases  of  interest  this  ratio  of 
element  diameters  becomes  unbounded  on  the  family  of  finite  element 
spaces  being  used.   In  this  case,  while  results  for  quasi-uniform 
spaces,  such  as  those  in  the  last  chapter,  yield  an  0(N)  work  bound  for 
each  finite  element  space,  the  constant  hidden  in  the  0(N)  work  bound 
will  depend  on  the  ratio  of  element  diameters.   Such  bounds  are  clearly 
unsatisfactory  for  locally  refined  finite  element  spaces. 

While  nearly  all  of  the  convergence  results  in  the  literature 
apply  only  to  the  case  of  quasi-uniform  finite  difference  or  finite 
element  grids,  multi-level  methods  are  being  used  increasingly  for 
locally  refined  grids,  Brandt  (1977),  Bank  and  Sherman  (1978).   In  fact, 
one  of  their  principal  advantages  is  that  their  rate  of  convergence 
remains  very  good  on  locally  refined  grids.   By  contrast  other  iterative 
methods,  whose  rates  of  convergence  depend  on  the  condition  number  of 
the  stiffness  matrix,  do  poorly  on  locally  refined  grids.   For  example, 
the  number  of  iterations  required  with  conjugate  gradient  iteration 


goes  as  /K(K),  where  K(K)  is  the  spectral  condition  number  of  the  stiffness 
matrix,  Axelsson  (1977).   This  condition  number  satisfies 

K(K)  >  c  h"2 
—    mm 

where  h  .   is  proportional  to  the  diameter  of  the  smallest  element  in  the 
mm 

grid.   Thus,  K(K)  may  grow  much  faster  as  a  function  of  the  number  of 
unknowns  on  a  locally  refined  grid.   This  presents  a  very  severe 
difficulty  for  nearly  all  iterative  methods.   As  will  be  shown  in  this 
chapter,  multi-level  methods  do  not  suffer  from  this  difficulty. 

The  algorithms  in  this  chapter  differ  considerably  from  that  in 
the  last  chapter,  aside  from  their  application  to  locally  refined  grids. 
In  the  algorithm  of  the  last  chapter,  to  approximately  solve  the  elliptic 


38 

problem  on  the  finest  grid,  a  related  problem  on  the  next  coarser  grid 
was  solved  twice.   It  was  solved  once  at  the  beginning  to  obtain  a 
starting  value  for  the  simultaneous  displacement  iteration  and  once  at 
the  end  as  a  final  correction.   In  the  algorithms  of  this  chapter  and 
in  other  algorithms  in  the  literature,  solution  of  a  related  problem  on 
the  next  coarser  grid  is  used  more  frequently,  typically  after  every 
few  smoothing  iterations  on  the  finest  grid.   Thus,  the  basic  operation 
in  such  an  algorithm  is  best  thought  of  as  a  multi-level  iteration 
consisting  of  a  simple  relaxation  iteration  periodically  accelerated 
by  the  solution  of  a  related  problem  on  the  next  coarser  grid. 

Algorithms  like  those  considered  in  this  chapter  have  been 
described  for  quasi-uniform  grids  by  Bank  and  Dupont,  Nicolaides, 
Hackbusch,  and  others.   They  are  also  related  to  methods  for  locally 
refined  grids  suggested  by  Brandt  and  very  closely  to  the  method  used  in 
the  adaptive  finite  element  code  of  Bank  and  Sherman.   An  algorithm 
similar  to  that  of  chapter  three  could  have  been  used,  but  probably 
would  not  have  been  as  practical  as  the  algorithms  here.   One  could  argue 
that  in  any  case  it  is  best  to  consider  different  modifications  of  the 
standard  approach  separately  rather  than  in  combination. 

To  our  knowledge  there  is  only  one  convergence  result  in  the 
literature  for  multi-level  methods  on  locally  refined  grids.   This  is 
given  in  Bank  and  Dupont  (1978)  where  a  two-level  multi-level  iteration 
is  shown  to  converge  at  a  rate  independent  of  the  mesh  size.   Since  this 
result  relies  only  on  local  properties  of  the  finite  element  space  and 
not  on  global  approximation  properties,  it  applies  automatically  to 
locally  refined  grids.   The  limitation  to  this  result  is  that  it  deals 
only  with  their  two-level  scheme.   With  such  a  scheme,  as  one  decreases 
the  mesh  size,  the  dimensions  of  both  the  finer  and  coarser  finite 
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element  spaces  grow.   In  their  algorithm  the  linear  system  for  the 
coarser  finite  element  space  must  be  approximately  solved  once  per 
multi-level  iteration.   Whether  this  is  done  directly  or  by  one  of 
the  common  iterative  methods,  the  cost  of  approximately  solving  this 
linear  system  grows  faster  than  linearly  in  the  dimension  of  the 
linear  system.   As  a  consequence,  such  a  two-level  scheme  can  never  be  of 
optimal  order.   It  may,  however,  be  quite  attractive  as  a  practical 
algorithm  since  inversion  of  the  lower  dimensional  linear  system  for 
the  coarser  space  may  be  much  cheaper  than  the  cost  of  inverting  the 
linear  system  for  the  finer  space  in  which  the  solution  is  sought. 
The  use  of  direct  methods  to  solve  these  linear  systems  for  the  coarser 
space  is  particularly  attractive  since  after  the  first  multi-level 
iteration  only  back  solves  are  required. 

One  might  expect  that  a  simple  generalization  of  their 
convergence  proof  for  the  two-level  case  could  be  used  to  show  optimal 
order  convergence  of  a  multi-level  method  using  many  levels  on  locally 
refined  grids.   In  some  cases  it  can,  for  example,  for  Poisson's  equation 
on  certain  finite  element  grids.   In  general,  however,  the  error 
reduction  per  iteration,  y   E  (0,  1),  of  their  two-level  scheme,  while 
independent  of  the  mesh  size,  depends  strongly  on  the  PDE  and  on  certain 
properties  of  the  discretization.   The  inductive  argument  used  to  show 
the  optimal  order  convergence  of  multi-level  methods  using  many  levels 
requires  that  we  be  able  to  choose  a  sufficiently  small  value  of  y. 
For  example,  for  the  first  algorithm  in  this  chapter,  y   must  be  chosen 
in  (0,  — ) .   Since  there  is  no  easy  way  of  controlling  y   in  their  two-level 
scheme,  it  seems  to  be  quite  hard  to  get  general  optimal  order  convergence 
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results  from  a  direct  generalization  of  their  two-level  result.   This  is 
unfortunate,  since  the  argument  showing  the  convergence  of  their 
two-level  algorithm  is  far  simpler  than  that  given  here. 

The  results  here  and  the  technique  of  their  proof  are  actually 
much  more  closely  related  to  the  second  kind  of  algorithm  considered  in 
their  paper.   This  second  algorithm  is  a  recursive  algorithm  using  many 
grid  levels  just  as  the  algorithm  in  the  last  chapter  did.   Their 
algorithm  is  shown  to  converge  in  0(N)  operations  on  quasi-uniform 
finite  element  spaces.   Unlike  the  results  in  chapter  three,  their 

results  apply  also  to  non-convex  domains,  on  which  the  elliptic  problem 

2 
is  less  than  H  regular.   The  analysis  of  algorithms  for  locally  refined 

grids  here,  which  also  applies  to  non-convex  domains,  is  largely 

modeled  on  their  analysis  of  the  second  algorithm  in  their  paper.   The 

notations  here  have  also  been  kept  almost  the  same  as  theirs. 

Analysis  of  multi-level  methods  is  harder  than  that  of  most 

iterative  methods  since  it  involves  both  algebraic  and  approximation 

questions.   The  algebraic  and  complexity  questions  for  locally  refined 

grids  are  slightly  harder,  since  the  dimensions  of  the  finite  element 

spaces  used  in  the  multi-level  algorithm  must  be  taken  into  account  more 

carefully.   By  contrast,  the  approximation  questions  for  locally  refined 

grids  seem  to  be  a  great  deal  harder.   In  the  first  half  of  this  chapter 

we  consider  the  algebraic  and  complexity  questions  involved  in  treating 

multi-level  methods  on  locally  refined  grids.   The  analysis  here 

largely  follows  that  in  the  paper  by  Bank  and  Dupont,  while  the  algorithms 

themselves  are  motivated  by  the  adaptive  finite  element  code  of  Bank 

and  Sherman.   Two  algorithms  are  considered;  the  first  is  similar  to 

that  in  the  Bank  and  Sherman  code.   The  second  is  a  modification  yielding 

an  improved  complexity  bound.   This  modification,  which  appears  to  be 
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new,  allows  greater  freedom  in  the  choice  of  finite  element  spaces 
used  and  may  have  significant  practical  applications. 

The  second  half  of  this  chapter  deals  with  the  approximation 
questions  involved  in  treating  locally  refined  grids.   The  approximation 
result  proved  there  is  the  main  result  of  this  thesis.   In  previous  work 
on  multi-level  methods,  the  approximation  questions  were  answered  by 
using  the  elliptic  regularity  of  the  problem  to  show  smoothness  of  the 
functions  being  approximated.   Then  the  approximation  error  could  be 
estimated  using  standard  finite  element  techniques  including  duality 
arguments  and  the  Bramble-Hilbert  lemma.   The  approximation  result  here 
uses  a  completely  different  kind  of  argument  making  no  use  of  the 

regularity  of  the  problem.   This  is  why  it  applies  to  non-convex 

2 
domains  where  the  regularity  is  weaker  than  H  .   Perhaps  of  greater 

interest  from  a  computational  point  of  view  is  the  fact  that  this  approx- 
imation result  relies  only  on  local  properties  of  the  finite  element 
spaces  and  the  PDE  operator.   In  consequence, it  would  be  quite  easy  to 
compute  rigorous  upper  bounds  on  the  rate  of  convergence  of  the  algorithms 
in  this  chapter.   Up  until  now  there  were  good  heuristic  bounds  on  the 
rate  of  convergence  of  multi-level  methods  on  finite  difference  grids, 
Brandt  (1977),  and  in  a  few  cases  rigorous  bounds.   No  previous  bounds 
for  the  rate  of  convergence  of  multi-level  methods  on  irregular  finite 
element  grids  appear  to  be  available. 

4.2.   Notation 

To  develop  the  algorithms  considered  in  this  chapter,  it  is 
necessary  to  go  into  greater  detail  concerning  the  finite  element  spaces 
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used  than  was  required  in  the  last  chapter.   Additional  notation  must 
also  be  developed.   These  are  the  tasks  of  this  section.   To  keep  the 
mathematical  issues  in  this  chapter  clear,  our  analysis  will  be 
restricted  to  C   triangular  and  rectangular,  elements.   Though  the 
analysis  in  this  chapter  seems  to  be  extendable  to  all  or  nearly  all 
families  of  finite  elements,  the  details  of  the  proofs  here  can  be 
considerably  harder  for  some  families  of  elements.   Using  the  C   square 
and  triangular  elements  considered  here  permits  considerable  simplifi- 
cation of  the  analysis  and  may  also  cover  the  cases  of  greatest  interest, 
since  the  majority  of  finite  element  computation  seems  to  be  done  with 
these  elements.   Linear  C  triangular  elements,  called  Turner  triangles, 
are  used  in  the  Bank  and  Sherman  adaptive  multi-level  code,  so  at  least 
for  polygonal  domains,  the  solution  space  used  by  their  code  is  covered 
by  the  theory  here.   However,  as  will  be  seen,  their  algorithm  violates 
the  assumptions  of  the  theory  developed  here  in  a  number  of  other  ways. 
There  are  really  two  reasons  why  more  detailed  specification 
of  the  finite  element  spaces  used  in  the  multi-level  algorithm  is  needed 
here  than  was  required  in  the  last  chapter.   The  first  is  that  locally 
refined  spaces  are  intrinsically  more  complex  than  quasi-uniform  spaces 
since  they  are  not  characterized  by  a  single  mesh  parameter.   The 
second  is  that,  unlike  the  analysis  in  the  last  chapter,  which  was 
based  on  standard  finite  element  error  estimates,  the  analysis  here  is 
based  on  an  approximation  result  developed  here.   The  proof  of  this 
result  depends  on  detailed  analysis  of  the  finite  element  spaces  used. 
While  this  result  could  have  been  proven  more  abstractly,  the  approach 
here  is  certainly  more  straightforward.   The  reader  with  a  background 
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in  finite  element  theory  should  have  little  trouble  in  seeing  the 
generalizations  to  elements  with  higher  order  continuity. 

The  two  elements  considered  here  are  the  standard  C  triangular 
and  square  elements.   The  C  triangular  elements  of  degree  k,  k  _>  1, 
use  as  element  trial  space  the  space  P  of  polynomials  in  x  and  y 

K. 

of  degree  at  most  k.   The  C  square  elements  of  degree  k,  k  >_  1,  are 
tensor  product  elements  based  on  the  trial  space  (L  consisting  of  poly- 
nomials  whose  terms  are  of  degree  k  or  less  in  x  and  y  separately.   Both 
types  of  elements  may  be  taken  as  nodal  finite  elements  with  node 
placement  as  in  figures  1  and  2. 
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Figure  1.   C  Lagrangian  Triangular  Elements 
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Figure  2.   C  Tensor  Product  Elements 
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The  nodal  parameters  used  are  just  the  function  values  at  these  points, 
providing  a  substantial  simplification  of  the  analysis  in  this  chapter 
over  that  which  would  be  required  for  the  more  complex  Hermite  elements. 
Equally  important  from  our  point  of  view,  for  any  k  _>  1 ,  either  type 
of  element  constitutes  an  affine  family,  Ciarlet  (1979).   What  this  means 
is  that  we  may  fix  a  reference  element  in  the  plane,  on  which  we  have 
a  reference  trial  space  L,  which  will  be  either  V     or  ()   for  triangular 

K.       K. 

or  square  elements,  respectively.   For  any  other  element  of  the  same 
type  anywhere  in  the  plane  we  can  find  an  affine  transformation  from 
the  reference  element  to  the  given  element.   This  affine  transformation 
carries  the  region  occupied  by  the  given  element,  the  nodes,  and  the 
element  trial  space  from  the  reference  element  to  any  other  element 
of  the  same  type  anywhere  in  the  plane. 

Affine  transformations  are  mappings  of  the  form 

x1  =  Ax  +  b 
where  b  is  a  fixed  coordinate  vector  and  A  is  a  nonsingular  2><2 
matrix.   For  triangular  elements  all  affine  transformations  may  be 
used,  carrying  the  reference  element  into  the  entire  six  parameter 
family  of  triangles  in  the  plane.   For  square  elements,  A  may  be  taken 
as  a  multiple  of  the  identity.   That  is,  there  is  no  need  for  rotation 
and  distortion  in  this  case. 

The  point  of  all  this  is  that  affine  transformations  and  the 
use  of  reference  elements  provide  an  elegant  simplification  of  the  theory. 
The  reference  element  and  its  element  trial  space  basically  contain  all 
properties  intrinsic  to  this  type  of  finite  element.   On  the  other 
hand,  the  affine  transformations  contain  all  information  about  the 
geometry,  mesh  size,  and  distortion  of  any  particular  element  in  the 
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finite  element  grid.   Thus,  analysis  of  a  finite  element  space  is 
conveniently  factored  into  two  parts,  a  property  we  will  use  to 
advantage  in  section  4.5. 

Aside  from  the  choice  of  element,  two  types  of  grids  need  to 
be  considered.   These  may  be  conveniently  labeled  "aligned"  and 
"unaligned"  grids.   Aligned  grids  are  those  which  are  standard  in  finite 
element  theory,  where  each  edge  of  an  element  is  the  edge  of  another, 
or  part  of  the  boundary.   Unaligned  grids  are  those  where  the  edge  of  an 
element  may  be  only  part  of  the  edge  of  another.   See  figures  3,  4,  and  5 
for  examples  of  these  two  types  of  grids.   In  general,  triangular  grids 
may  always  be  chosen  to  be  aligned,  without  sacrificing  efficiency,  but 
it  is  often  necessary  to  choose  grids  of  square  or  rectangular  elements 
to  be  unaligned,  or  much  of  the  benefit  of  local  refinement  will  be  lost. 

To  be  admissible,  the  finite  element  spaces  must  be  C  and 
satisfy  the  essential,  that  is  Dirichlet,  boundary  conditions.   On 
aligned  grids  there  is  no  problem  obtaining  C  continuity.   One  simply 
requires  agreement  of  the  nodal  parameter  on  the  edge  between  adjacent 
elements.   This  does  not  work  for  unaligned  grids.   To  simplify  this 
problem  we  restrict  the  unaligned  grids  as  follows.   Assume  that  every 
edge  of  an  element  is  either  part  of  the  boundary,  the  edge  of  another 
element,  or  the  union  of  the  edges  of  two  other  elements.   This  last 
situation  is  shown  in  figure  6. 

In  this  last  situation  we  consider  the  single  element  to 
dominate  the  two  elements  it  adjoins  across  this  edge.   For  convenience 
call  the  single  element  the  "larger"  since  it  almost  always  is  larger 
than  the  other  two.   Then  we  may  say  that  the  "larger"  dominates  the 


46 


Figure  3.   Aligned  Triangular  Grid 
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Figure  4.   Unaligned  Triangular  Grid 


Figure  5.   Unaligned  Grid  with  Square  Elements 
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Figure  6.   Adjacent,  Unaligned  Elements 


47 

other  two  in  the  sense  that  its  nodes  on  this  edge  are  the  real  nodes 

of  the  finite  element  space  there.   The  nodes  of  the  "smaller"  elements 

may  be  called  parasitic  since  their  values  will  be  taken  from  the  larger 

element.   Parasitic  nodes  of  the  "smaller"  elements  in  figure  6  not 

coinciding  with  nodes  of  the  "larger"  element  are  labeled  with  a  "P." 

It  is  easy  to  see  that  this  system  preserves  the  required  C  continuity 

and  is  fairly  straightforward  to  code. 

One  more  item  needs  to  be  discussed  before  we  consider 

construction  of  the  locally  refined  spaces  needed  in  this  chapter.   This 

is  the  angle  condition.   It  was  mentioned  in  section  2.3  that  the  mass 

matrix  will  usually  be  well  conditioned  uniformly  in  the  level  of 

refinement  of  the  grid  as  long  as  the  element  geometry  does  not 

degenerate.   While  this  is  not  always  true  for  locally  refined  grids, 

it  is  still  important  that  the  element  geometry  not  degenerate.   One 

way  of  stating  the  required  condition  is  to  say  that  the  minimum  interior 

angle  of  any  triangle  in  the  grid  must  be  bounded  below  uniformly  in 

the  level  of  refinement.   More  formally,  if  {T.j.  ,  is  the  family  of 

3  3=1 

triangulations  being  considered,  we  want 

(4.2.1)  inf  {  min  <J>  }  >  0 

j>l  T^T. 
J 

where  6   is  the  smallest  interior  angle  of  the  triangle  T.   There  are 
several  alternate  formulations  of  this  condition  in  the  finite  element 
literature,  Oden  and  Reddy  (1976),  but  this  one  is  adequate  for  our 
purposes.   Obviously,  no  such  condition  is  needed  on  square  grids. 

We  next  turn  to  consideration  of  the  class  of  finite  element 
grids  on  which  the  multi-level  algorithms  of  this  chapter  can  be  shown 
to  converge  rapidly.   This  class,  which  includes  both  locally  refined  and 
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quasi-uniform  spaces,  is  quite  large  and  may  be  thought  of  as  an  indexed 

family  {M  }  _   of  spaces.   It  consists  of  all  spaces  which  can  be 
a  oea 

generated  by  the  refinement  process  described  below  and  which  satisfy 
an  additional  condition  given  in  section  4.4.   An  0(N)  complexity  bound 
holds  uniformly  on  the  class  of  spaces  in  {M  )__.  satisfying  this 

Ot  OfcA 

condition  uniformly.   In  general,  this  class  of  spaces  with  a  uniform 

0(N)  complexity  bound  is  large  enough  to  show  the  improved  convergence 

to  singular  solutions  obtainable  by  discretization  with  locally  refined 

grids.   A  partial  exception  to  this  is  the  case  of  point  singularities, 

since  the  best  discretizations  for  these  problems  do  not  satisfy  the 

conditions  of  section  4.4  uniformly.   For  these  problems  0(N)  bounds 

can  be  shown  by  the  theory  here  only  for  suboptimal,  but  still  very  good 

locally  refined  grids. 

Fixing  a  coarse  space  M  ,  the  refinement  process  described 

below  will  generate  sequences  {M  „  . } .  , ,  3  e  B  of  finite  element  spaces, 

3,1  3=1 

where  the  index  3  ranges  over  the  uncountable  family  of  choices  allowed 
in  this  refinement  process.   The  family  of  spaces  {M  }    is  just  the 
union  of  all  spaces  which  can  be  produced  by  this  refinement  process, 

The  condition  mentioned  above,  which  is  described  more  fully  in  section  4.4, 

is  that  for  each  fixed  3  G  B,  the  dimensions  of  the  spaces  {M0  ./.  . 

3,3  3=1 

in  this  sequence  of  refinements  grow  geometrically.   The  0(N)  complexity 
bound  given  here  depends  only  on  the  problem,  the  coarsest  finite 
element  space  M  ,  and  this  geometric  growth  rate. 
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The  multi-level  algorithms  here  actually  use  only  the  finite 

family  of  spaces  {M0  .}.  ,,  for  fixed  3  e  B  and  fixed  k  >  2,  in  their 

3,j  3=1'  - 

computation.   Since  3  is  fixed,  it  is  convenient  to  drop  it,  although  we 

oo 

retain  j,  and  in  fact  consider  the  entire  sequence  of  refinements  (M.}._ 

since  that  simplifies  the  statements  of  some  of  the  theorems  here.   One 

word  of  warning  is  needed,  however.   Though  the  spaces  in  the  sequence 

{M.}.  n    become  infinite  dimensional,  they  do  not  usually  become  dense  in 

H'  and  thus  convergence  to  the  true  solution  will  not  ordinarily  be 
£ 

oo 

obtained  in  {M.}._  .   In  order  to  get  a  sequence  of  solutions  converging 
to  the  true  solution,  one  must  go  back  to  consideration  of  the  family 

a  a£A 

OO 

Description  of  the  families  (M . }    of  locally  refined  grids 
is  essentially  the  same  for  grids  based  on  square  and  triangular  elements. 
We  consider  first  the  triangular  case  since  it  is  the  case  of  greatest 

interest  in  this  chapter.   In  this  case,  the  family  {M.}._,  of  finite 

r   i00 
element  spaces  will  be  based  on  a  nested  family  of  triangulations  lT.j.  ... 

By  nested  we  mean  that  each  triangle  T  of  T.,  j  >  2,  is  contained  in  a 

r   ">  °° 

triangle  of  T.  ,.   Once  we  have  chosen  the  triangulations  It.;.  ..  ,  and 
3-1  J  J=1 

the  particular  element  to  be  used,  the  family  of  finite  element  spaces 

r   i  °° 

W.i.    .    is  completely  determined. 

J  3=1 

The  construction  of  the  nested  family  of  triangulations  It  j.  , 

j  J=1 

may  be  done  by  a  process  Bank  and  Sherman  call  regular  subdivision.   In 

regular  subdivision  of  a  triangle  T,  it  is  divided  into  four  smaller  triangles 

by  pairwise  connecting  the  midpoints  of  the  edges  together  as  shown  in 

figure  7.   Let  T  be  a  fixed,  aligned  triangulation  of  Si.      Now  inductively 

suppose  T.  has  been  chosen  for  some  i  >  1.   Select  a  subset  T,  of  x..   We 
3  ~  J     J 
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may  take  T    as  the  (unaligned)  triangulation  of  fi  formed  by  regular 
subdivision  of  the  triangles  in  T.  while  leaving  the  remaining  triangles 
of  T  unaltered.   Proceeding  is  this  way,  the  entire  family  of 


oo 

triangulations  {t,}.  ,  may  be  constructed, 
J  3=1 


Figure  7.   Regular  Subdivision  of  a  Triangular  Element 

In  our  approach  the  subsets  T.  C  t  ,  j  >  1,  of  triangles 

to  be  refined  cannot  be  chosen  arbitrarily,  but  most  satisfy  certain 

conditions  necessary  for  the  operation  and  analysis  of  the  algorithms 

here.   Let  £2.,,  C  fi  be  the  subdomain  convered  by  triangles  in  T.,  for 
j+l  J 

j  _>  1.   For  simplicity  of  notation  it  is  convenient  to  set  (L  =  fi. 

The  first  condition  required  is 

(4.2.2)  0    CD,  ,  j  >  1  . 

J+l    J     - 

This  condition  means  that  for  j  >  1,  only  triangles  of  T.  which  were 
formed  by  regular  subdivision  of  triangles  in  T.  -  may  be  subdivided 
in  creating  T .  _ .  Equivalently,  for  j  _>  1»  each  triangle  of  T.  -  is 
contained  in  a  triangle  of  T,. 

Now  assume  that  the  diameters  d™  of  the  triangles  T  in  the 
coarsest  triangulation  T   satisfy 

(A. 2. 3)  -  hn  <  d  <  a  h,  , 

a  1  —  T  —    1 

for  some  fixed  a  >  1  and  mesh  parameter  h  .   As  before  set 
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(4.2.4)  h.  =  p1'2    h±    ,  j  >  1  , 

where  p  is  a  mesh  ratio,  which  for  the  regular  subdivision  process  here 
is  two.   Conditions  (4. 2. 2)- (4. 2. 4)  imply  that  the  diameters  of 
triangles  in  T .  ,  j  >_  1,  must  satisfy 

(4.2.5a)  -  h.  <  d   <  a  h.  for  T  C  ft. 

o     J  -  T  -    j  j 

(4.2.5b)  ~  h.  <  d_  <  a)  h.  for  T  C  ft.  \  n  ., 

a  i  —  T  —    i         i   i+I 

with  1  <  i  <  j  -  1  . 
The  overbar  here  denotes  closure,  necessary  since  triangles  are  closed, 
and  our  domains  {ft.}.  ,  are  not. 

Intuitively,  the  triangles  in  T  are  all  of  about  the  same 
size.   By  (4.2.3)  or  (4.2.5a)  their  diameters  are  all  comparable  to  h  . 
The  triangles  in  T   come  in  two  sizes:   fine  ones  in  ft~  with  diameter 
comparable  to  h_,  produced  by  regular  subdivisions,  and  coarse  ones 
outside  ft„  with  diameters  about  h  .   In  general,  T.  may  be  thought  of 
as  having  j  sizes  of  triangles,  the  finest  lying  in  ft..   Of  course, 
this  intuitive  point  of  view  should  not  be  pushed  too  far.   Because  of 
the  constant  a  in  (4.2.5)  these  sizes  may  overlap  a  great  deal. 

Only  one  other  condition  is  really  essential  for  triangular 
grids.   That  is,  that  the  number  of  elements  adjoining  a  given  element 
along  one  edge  be  at  most  two.   The  "two"  here  is  arbitrary,  and  any 
other  finite  bound  could  be  chosen,  but  in  practice  two  is  probably  all 
one  wants.   The  coding  is  difficult  enough  as  it  is.   This  restriction, 
that  at  most  two  "smaller"  elements  adjoin  one  "larger"  element  across 
an  edge,  together  with  the  angle  condition  (4.2.1)  limits  the  rate  of 
change  of  the  mesh  size  as  one  can  easily  show. 
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Locally  refined  grids  of  square  elements  may  be  constructed  along 
the  same  lines  as  the  triangular  grids  just  considered.   One  begins  with 

an  aligned  square  grid  G     on  ft  and  then  generates  a  family  of  finer 

.00 
(unaligned)  grids  (G.}    ,  using  an  analogous  regular  subdivision  process. 

This  process  is  shown  in  figure  8.   The  subdomains  (ft./.  ,  can  be  defined 

analogously,  and  the  analogs  of  (4. 2. 3)-(4. 2. 5)  will  hold,  now  with  the 

notation  "d  "  instead  of  "d  ",  using  "e"  for  "element."   Since  such 
e  T       ° 

square  grids  are  only  being  considered  in  this  section  and  will  thereafter 
be  dropped  for  simplicity,  there  is  no  need  to  go  into  these  details 
more  deeply. 


Figure  8.   Regular  Subdivision  of  a  Square  Element 
A  problem  arises  on  locally  refined  grids  not  present  on 
quasi-uniform  grids.   Unless  special  care  is  taken,  the  condition  numbers 
of  the  mass  matrices  for  a  family  of  locally  refined  finite  element  spaces 
may  not  be  uniformly  bounded.   For  aligned  quasi-uniform  grids  and  the 
usual  types  of  elements,  including  those  considered  here,  the  uniform 
boundedness  of  the  condition  number  of  the  mass  matrices  is  equivalent 
to  the  angle  condition  (4.2.1).   Unaligned  quasi-uniform  grids  are  rarely 
considered,  for  obvious  reasons,  but  the  same  thing  can  be  shown  for 
such  grids. 
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The  problem  with  locally  refined  grids  is  that  the  number  of 
basis  functions  which  overlap  with  a  given  basis  function  may  be  very 
large,  becoming  unbounded  as  the  grid  is  refined.   This  violates  the 
hypotheses  of  the  result  of  Fried  described  in  section  2.3,  and  in  fact 
mass  matrices  with  condition  numbers  becoming  unbounded  do  occur.   A 
simple  finite  element  grid  where  this  problem  occurs  is  shown  in 
figure  9.   Here  bilinear  elements  are  used  with  the  standard  "hat  function" 
or  "pagoda"  basis  functions.   Notice  that  the  conditioning  problem  in 
this  example  occurs  despite  the  fact  that  the  mesh  size  changes  slowly 
in  the  sense  that  the  diameters  of  adjacent  elements  differ  at  most  by 
a  factor  of  two.   Proper  scaling  of  the  basis  functions  mitigates  the 
problem  but  does  not  eliminate  it. 


sei; 

i-+ 

++ 

f+ 

i 

f 

f+ 

+i 

+  4 

Figure  9.   Finite  Element  Space  with  Poorly  Conditioned  Mass  Matrix 

The  only  real  solution  is  to  arrange  things  so  that  each  basis 
function  overlaps  only  a  bounded  number  of  elements,  with  the  bound 
independent  of  the  level  of  refinement.   This  can  be  done  either  by 
modifying  the  basis  elements  to  prevent  excessive  overlap,  or  by  selecting 
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only  grids  where  such  overlap  does  not  occur.   The  second  approach  is 
probably  preferable  in  most  cases  since  it  tends  to  be  simpler.   Consider 
the  approach  of  modifying  the  basis  functions  for  the  simple  bilinear 
elements.   The  situation  is  similar  for  the-  other  square  and  triangular 
elements  here.   Ordinarily,  the  nodal  parameter  at  the  intersection  of  the 
four  elements  shown  in  figure  10  would  be  associated  with  the  "hat 
function"  whose  support  is  the  union  of  the  four  elements.   When  one  or 
more  of  these  four  elements  is  refined,  the  basis  function  associated  with 
its  center  node  can  be  modified  to  have  support  on  a  smaller  region. 
Up  to  rotation,  there  are  four  cases  to  consider,  as  shown  in  Figure  11. 
Clearly  the  case  where  all  four  elements  are  refined  need  not  be  dealt 
with  since  then  the  same  considerations  apply  on  the  resulting  smaller 
elements.   The  dark  lines  in  figure  11  show  the  boundaries  of  the 
supports  of  the  modified  basis  functions.   The  squares  marked  with 
asterisks  may  be  subject  to  further  refinement  without  violating  the 
requirement  that  the  mesh  size  change  slowly.   The  point  is  that  all 
squares  so  marked  lie  outside  the  support  of  these  basis  functions,  so 
the  number  of  elements  these  basis  functions  overlap,  and  hence  the 
number  of  other  basis  functions  they  overlap,  will  be  bounded  by  a  fixed 
constant  independent  of  the  level  of  refinement.   The  only  difficulty  with 
this  approach  is  that  computations  of  the  element  integrals  is  somewhat 
more  complicated  since  these  basis  functions  are  more  complex.   This 
approach  was  suggested  to  the  author  by  D.  Gannon,  who  has  done  considerable 
work  with  rectangular  tensor  product  elements  on  locally  refined  grids, 
Gannon  (1980). 


55 


<» 


Figure  10.   Support  of  a  Basis  Function 
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Figure  11.   Modified  Basis  Functions 
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In  the  second  approach,  one  retains  the  simple  Lagrangian  basis 
functions,  "hat  functions,"  and  their  higher  order  analogues,  but 
restricts  or  modifies  the  grids  so  that  the  problem  of  one  basis  function 
overlapping  an  increasingly  large  number  of  others  never  arises.   With 
triangular  grids  modification  is  not  difficult.   This  is  the  approach 
used  in  the  Bank  and  Sherman  code.   Beginning  with  an  unaligned 
triangulation  T . ,  1  _<  j  <  °°,  belonging  to  a  nested  family  of  triangulations 
IT.;.  ,  satisfying  the  conditions  given  so  far,  their  algorithm  constructs 
an  aligned  triangulation  x.,  refining  T.,  on  which  this  condition  number 
problem  does  not  arise.   This  is  done  in  two  steps.   In  the  first  step 
each  triangle  with  the  property  that  there  are  vertices  of  other 
triangles  in  the  middle  of  more  than  one  of  its  edges  is  regularly 
refined,  as  in  figure  12.   Once  all  such  triangles  are  refined,  including 
those  produced  in  carrying  out  this  process,  the  only  triangles  in  the 
grid  are  those  with  either  one  vertex  in  the  middle  of  an  edge  or  no 
vertices  in  the  middle  of  its  edges.   The  former  are  refined  by  a 
different  subdivision  process,  called  "green"  refinement,  shown  in 
figure  13.   The  resulting  aligned  triangulation  T.  is  then  used  to  generate 
the  finite  element  space  M.  used  by  the  multi-level  algorithm. 

Notice  that  green  refinement  may  decrease  the  smallest  interior 
angle  present  on  the  grid.   This  would  be  a  problem  if  successive  green 

refinements  were  applied.   In  Bank  and  Sherman's  approach  this  is  avoided 

r   i  °° 
by  first  generating  the  unaligned  triangulations  Ix.j    and  then 

r  ~  1°° 
forming  the  corresponding  aligned  triangulation  It  .  J- . _  ,  rather  than 

taking  T    ,  j  >  1,  as  a  refinement  of  x..   This  ensures  that  the  angle 

condition  (4.2.1)  is  satisfied.   Even  so,  the  triangulations  IT./.*  used 

by  their  code  are  not  generally  admissible  in  the  theory  here,  even  when 
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A 


Figure  12.   Regular  Refinement  of  Triangles  with  Edges  Containing  Vertices 


Figure  13.   Green  Refinement  of  a  Triangle 
the  underlying  unaligned  triangulations  It./.,  satisfy  all  of  our 
assumptions.   The  problem  is  that  in  producing  an  aligned  triangulation 
T .  ,  j  >_  1,  a  triangle  may  be  green  refined  while  in  a  finer  unaligned 
triangulation  T.,  i  >  j,  this  same  triangle  may  be  regularly  refined. 

r~  1°° 

If  so,  the  triangulations  It./    will  not  be  nested,  nor  will  the  resulting 

r   1°° 

finite  element  spaces  {M.i.    ,.   While  this  does  not  seem  to  be  a  problem 

J  J=l 

for  their  code,  it  is  certainly  a  problem  for  the  theory  here.   The  cure 

for  this  seems  to  be  to  require  the  unaligned  triangulations  It.} 

to  contain  only  triangles  with  a  vertex  of  another  triangle  in  at  most 
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one  of  their  edges.   If  so,  in  going  from  {x  }    to  the  corresponding 

r  ~  -.00 

aligned  triangulations  It)    only  green  refinement  will  be  done.   Since 


the  triangles  of  T  ,  j  >  1,  green  refined  in  producing  T.  must  lie 
outside  fi.  (the  portion  of  T.  lying  in  fi.  is  aligned),  there  is  no 
danger  of  subsequently  regularly  refining  them.   Further,  it  is  easy 
to  see  that  those  triangles  of  T.  produced  by  green  refinement  must  be 

green  refined  in  producing  T  ,  i  >  j ,  so  the  finite  element  spaces  based 

r~-|°° 
on  the  family  of  triangulations  it).,  will  be  nested  and  the  theory  of 

this  chapter  will  go  through  for  them. 

The  advantage  of  this  approach  is  that  the  condition  number 

of  the  mass  matrix  for  a  finite  element  space  based  on  an  aligned 

triangulation  will  be  uniformly  bounded,  as  long  as  the  angle  condition 

is  satisfied.   This  is  so  for  linear  elements  since  each  pyramid  basis 

function  overlaps  only  those  elements  meeting  at  the  vertex  with  which 

its  nodal  parameter  is  associated.   The  support  is  thus  a  polygon  formed 

from  these  elements,  as  in  figure  14.   The  other  kinds  of  basis  functions 

required  for  higher  order  elements,  those  associated  with  nodes  on  the 

edges  or  interiors  of  triangles,  are  not  a  problem  either,  since  their 

supports  are  smaller  than  those  associated  with  the  vertices. 


Figure  14.   Support  of  a  Basis  Function  on  an  Aligned  Triangulation 
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Though  the  above  approach  is  quite  attractive  for  triangular 

grids,  since  aligned  triangulations  are  convenient  in  programming,  it 

does  not  generalize  directly  to  square  elements,  since  there  is  no 

analog  of  green  refinement  for  these  elements.   Instead,  we  must  either 

use  the  modified  basis  functions  described  above,  or  restrict  the  rate 

of  change  of  the  mesh  size  so  that  the  conditioning  problem  does  not 

occur  when  using  the  standard  basis  functions.   The  simplest  way  to 

describe  the  required  restriction  on  the  rate  of  change  of  the  mesh 

size  is  in  terms  of  the  process  for  constructing  the  grids  iG .  J- ._  of 

square  elements  we  wish  to  consider.   Say  that  two  elements,  e  ,  e„  £  G., 

1  ■<  j  <  oo }  touch  if  they  have  any  point  in  common;  that  is,  if  they 

share  a  vertex  or  part  of  an  edge  as  in  figure  15.   Let  G\  ,  j  _>  1,  be 

the  set  {e  G  G.:  e  C  ft.}.   That  is,  G\    is  the  set  consisting  of  the 
J       J  J 

most  refined  elements  of  G..   Define  the  interior  of  G\ ,  denoted 

3  3 

int(G|)  as  the  set  of  elements  of  G\    not  touching  any  elements  in 

G.\  G\ .   Then  a  sufficient  condition  for  the  mass  matrices  for  the  locally 
3        3 

r   i  °°  r  r>     1°° 

refined  spaces  {M.}._  ,  based  on  the  grids  \G.}._   ,  to  be  uniformly 

well  conditioned  is 

G.  C  int(G|)  ,  j  >  1  . 
j        J      ~ 

This  condition  suffices  for  the  mass  matrices  to  be  uniformly  bounded 

since  no  basis  functions  centered  at  nodes  outside  of  G\    have  supports 

3 

extending  into  int(G'. ).   It  follows  that  they  only  overlap  elements  outside 
ft.,  and  so  they  only  overlap  a  bounded  number  of  elements. 
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'        <     1 


Figure  15.   Touching  Elements 
With  the  two  types  of  elements  considered  here  and  the 
variety  of  ways  of  ensuring  the  uniform  boundedness  of  the  condition 
numbers  of  the  mass  matrices,  a  number  of  options  are  available. 
However,  it  would  be  unwieldy  to  treat  the  theory  in  the  rest  of  this 
chapter  in  the  same  generality,  particularly  the  work  in  section  4.5, 
which  relies  on  detailed  analysis  of  the  finite  element  spaces  used. 
For  this  reason,  we  restrict  the  analysis  that  follows  to  C  spaces 
on  aligned  triangular  grids.   This  choice  is  motivated  partly  by 
convenience,  and  partly  by  the  use  of  such  spaces  in  the  Bank  and 
Sherman  code.   There  would  be  no  difficulty  in  treating  the  other 
types  of  finite  element  spaces  considered  here,  however,  the  interpo- 
lation mappings  of  section  4.5  are  slightly  more  complicated  to 
analyze  on  unaligned  grids. 
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4.3.   Algorithms 

In  this  section  we  describe  the  multi-level  algorithms  to  be 
considered  in  this  chapter  and  begin  the  analysis  of  their  convergence. 
The  first  of  these  algorithms  is  identical  to  one  of  those  described  by 
Bank  and  Dupont  for  quasi-uniform  grids  and  is  similar  to  algorithms 
described  by  Nicolaides,  Hackbusch,  Brandt,  Bokhvalov,  and  others.   It 
is  also  similar  to  the  algorithm  used  in  the  adaptive  finite  element 
code  of  Bank  and  Sherman  though  there  is  a  significant  difference  here 
which  will  be  discussed.   The  second  algorithm  to  be  considered  here 
appears  to  be  new.   It  can  be  used  to  show  the  0(N)  complexity  bound 
given  in  the  next  section  under  weaker  assumptions  than  otherwise  necessary, 

Both  of  the  algorithms  described  here  differ  substantially 
from  that  described  in  the  last  chapter.   There,  to  approximately  solve 


the  finite  element  equations  on  a  member  M.,  j  ^2,  of  a  family  {M.} 

J  ]  j  = 

of  finite  element  spaces,  approximate  solution  of  a  related  problem  in 
M 


_  was  required  twice;  once  at  the  beginning  to  obtain  a  starting 

value  for  smoothing  iteration  on  M.,  once  at  the  end  of  the  iteration  as 

a  final  correction.   The  algorithms  in  this  chapter  and  other  multi-level 

algorithms  in  the  literature  use  approximate  solution  of  the  equations 

for  M    more  often.   After  obtaining  a  starting  value  from  M._  ,  the 

accepted  solution  is  obtained  through  a  sequence  of  multi-level  (outer) 

iterations.   Each  of  these  outer  iterations  consists  of  a  number  of 

smoothing  (inner)  iterations  on  M.,  followed  by  approximate  solution  of 

a  related  problem  in  M.   .   Thus,  approximate  solution  of  the  related 

problem  in  M    is  required  an  indeterminate  number  of  times  depending 

on  how  far  the  outer  iteration  on  M.  is  carried. 

J 
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This  kind  of  algorithm  is  somewhat  more  complicated  than  that 
of  chapter  three.   However  it  made  sense  to  analyze  the  application  of 
this  kind  of  algorithm  to  locally  refined  grids  rather  than  the 
application  of  algorithms  like  those  in  the  last  chapter  for  several 
reasons.   The  main  one  was  the  desire  to  establish  the  optimal  order 
convergence  of  an  algorithm  as  close  as  possible  to  those  being  used 
in  practice.   This  was  more  fully  accomplished  than  it  would  have  been 
if  we  had  examined  algorithms  more  like  those  in  chapter  three.   It 
might  also  be  pointed  out  that  relatively  little  of  the  complexity 
of  this  chapter  is  caused  by  the  choice  of  algorithm. 

Description  of  the  algorithms  here  and  analysis  of  their 
convergence  properties  will  involve  the  eigenf unctions  and  eigenvalues. 
Following  the  notation  in  chapter  three  let  N.  be  the  dimensions  of  M., 

1  <.  j  <  °°.   Then  for  the  problem  (2.1.3)  we  have  eigenf  unctions 

(i)  Ni  (1)i^i 

U>)      } .    -,    in  M.  and  eigenvalues  {at/.  -  satisfying 
i   i=l     j  1  1=1 

aO^J),  <j0  =  A<j)(^j),  4>),  <J>  e  M.  . 

Since  more  than  one  finite  element  space  may  be  involved  in  our 
considerations,  we  will  retain  the  superscript  j.   As  in  chapter  three 
we  may  replace  the  L  inner  product  by  the  bilinear  form  b.(#,  •)  defined 
there,  where  we  now  use  the  subscript  j  to  keep  track  of  the  finite 


element  space  in  {M.}._   giving  rise  to  this  bilinear  form.   With  this 

-eplacemi 

{Xf^}.^    satisfying 
i        i=l 


replacement  we  have  as  before  eigenfunctions    {i|j.      } -Z.-,    and  eigenvalues 


a(^j),    cjO   =   A<j)    b.(^j),    ♦),   <Jre  Jl     . 

As  usual  we  assume  the  eigenvalues  are  numbered  in  order  of  increasing 
magnitude; 
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0  <  JL(J)  <  A<j)  <...<^}  , 

j 

o  <  3L(J)  <  A<j)  <...<  x<j)  . 

1   —  2   —   —  N. 

3 

For  both  of  the  algorithms  considered  in  this  chapter  the  basic 

multi-level  iteration  is  the  same.   In  both  cases  each  outer  iteration 

consists  of  a  sequence  of  smoothing  inner  iterations  followed  by 

approximate  solution  of  a  related  problem  in  a  coarser  finite  element 

space.   The  difference  between  the  two  algorithms  will  lie  in  the  manner 

of  solving  these  related  problems.   In  order  to  obtain  a  solution  in  a 

member  M  ,  k  >^  2,  of  a  family  (M.}.,  of  locally  refined  spaces, 

multi-level  iteration  will  need  to  be  applied  to  all  of  the  spaces 

{M.}.  .;  thus,  the  iteration  will  be  described  for  any  space  M.  in 
J  J=2  j 

3  3=1 

Beginning  with  a  trial  solution  u_   in  M.  the  outer  iteration 

0,n     j 

r  (i  )  V00 

produces  a  sequence  of  iterates  in   j    converging  to  the  true 

discrete  solution  u    in  M..   Each  of  these  outer  iterations  is 

J 

similar  to  the  algorithm  of  the  last  chapter.   Beginning  with  a  trial 

solution  u_   £  M.  an  outer  iteration  produces  an  improved  solution 
0,n    j 

U0^n+1  given  by 

u(j)   =  u(j) 
0,n+l    m+l,n  ' 

where  {unJ  }„  ,  is  a  sequence  of  inner  iterates.   All  but  the  last  of 
J6,n  ic=l 

these  inner  iterates  are  generated  by  a  simultaneous  displacement 

smoothing  iteration  on  M..   The  last  uses  a  correction  based  on  the 

3 

approximate  solution  of  a  related  problem  in  the  coarser  finite  element 

space  M.  . .   More  precisely,  beginning  with  a  trial  solution  un       £  M., 
j-1  0,n    j 

an  outer  iteration  consists  of  the  steps: 
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1.   For  I   =  1,...,  m  define  u„   £  M  by 


(a(u0  ,   ,  «J>)  -  (f,  W),46HJ( 


j 

2.   Letting  v  £  M._1  be  given  by 

(4.3.2)  a(v,  (j))  =  a(u   ,  $)  -  (f ,  <j>)  ,  <f>  €  U  , 

m,n  j-1 

let  v  be  an  approximation  to  v  and  set 

(4.3.3)  ii      =  u    -  v  . 

m+1 , n    m , n 

The  quality  required  in  the  approximation  of  v  by  v  will  be  specified 
in  theorems  4.1  and  4.3. 

The  intuitive  meaning  of  these  two  steps  is  that  in  the  first 
the  more  oscillatory  components  of  the  convergence  error  are  rapidly 
reduced  by  the  smoothing  iteration.   In  the  second  step  the  smoother 
components  of  the  error,  for  which  the  simultaneous  displacement 
iteration  is  ineffective,  are  reduced  through  approximate  solution  of  the 
related  problem  (4.3.2).   Consequently  convergence  will  be  rapid  for  all 
components  of  the  error. 

To  make  these  considerations  precise,  let  0 .  £  (0,  1)  be  an 
arbitrary  parameter.   Let  J (6.)  be  the  integer  such  that 

x(j)  <  e_  A(j)  for  i  <  j(e  ) 

i   —  3   N.        —    2 

3 

x(j)  >   e>  x(j)  for  ±  >  j(0j  . 
l      j   N. 

3 

Similarly  let  J(0.)  be  the  integer  such  that 
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A?j)  <  9.  X^j)  for  i  <  J(0.)  , 


i      N. 

L(j)  >  6.  \[p    for  i  >  J(6.)  . 
J 

Then  we  can  define  the  subspace  S.    0.   S.    V,    of  M.  by 

3      3      J   J     J 

S.  =  span  {\p^\    i  <  J(6.)}  , 

0.  =  span  ii>^\    i  >  J(6.)}  , 

S.  =  span  {$fj),  i  <  J(8.)}  , 
3  i     -    j 

0.   =   span  {ipfj),  i  >  J(6.)}  . 
J         i  J 

The  spaces  S.  and  S.  consist  of  relatively  smooth  functions  while  0. 
3  3  3 

and  0.    contain  more  oscillatory  functions.   For  any  9.  £  (0,  1)  one  has: 

M.  =  5.  ©  0.   =  S.    ©  0. 
3  3    ^      3  3    w   J 

In  fact,  0.    is  the  orthogonal  complement  of  5.  in  M.  in  both  the  energy 

and  L  inner  products,  while  0.    is  the  orthogonal  complement  of  S.  in 

3  3 

energy  and  the  inner  product  b^  (•,  •). 

For  the  time  being  only  the  subspaces  5.  and  V.   will  be  important 
in  analyzing  the  convergence  of  the  multi-level  iteration  (4. 3.1)-(4. 3.5) , 
since  the  approximation  question  involved  will  be  deferred  until  sections 
4.5-4.7.   There  S.  and  0.   will  play  an  important  role.   To  establish  the 
rapid  convergence  of  iteration  (4. 3.1)-(4. 3. 3)  we  will  show  that  error 
components  in  the  oscillatory  space  0.    are  rapidly  reduced  in  step  1, 
while  components  in  the  smooth  space  S.  are  reduced  in  step  2. 

The  multi-level  iteration  (4. 3.1)-(4. 3. 3)  differs  from  that 
in  the  Bank  and  Sherman  adaptive  finite  element  code  primarily  in  the 
smoothing  inner  iteration  used.   Their  code  uses  a  symmetric  Gauss-Seidel 
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smoothing  iteration  rather  than  the  simultaneous  displacement  iteration 
used  here.   While  comparisons  for  quasi-uniform  grids  have  shown  that 
the  choice  of  smoothing  iteration  is  not  very  important,  Brandt  (1977), 
Bank  and  Sherman  (1979),  for  locally  refined  grids  the  choice  is  significant, 
In  Fourier  terms,  the  iteration  (A.  3.1)  rapidly  annihilates  Fourier  error 
components  whose  wave  length  is  comparable  to  h  ,  the  finest  mesh  spacing 
of  M..   Away  from  ft.,  the  most  refined  region,  this  iteration  has  little 
effect.   The  reason  for  this  is  that  the  residual  on  the  right  hand 
side  of  4.3.1  is  multiplied  by  the  constant  factor  -  Z7T\    •      By  contrast 

y 

in  Jacobi  or  Gauss-Seidel  iteration  each  component  of  the  residual  is 
multiplied  by  the  inverse  of  the  corresponding  diagonal  entry  of  the 

stiffness  matrix.   In  other  words,  in  these  latter  cases  the  scaling 

2  2 

factor  looks  like  -ch,    ,  rather  than  -ch.,  where  h,    -  is  a  function 

local  j         local 

which  takes  on  a  different  value  for  each  component  of  the  residual 
vector,  depending  on  the  size  of  the  corresponding  elements.   As  a 
result  of  this  difference,  the  smoothing  iteration  used  in  the  Bank  and 
Sherman  code  is  effective  everywhere  on  the  grid,  rather  than  just  on  the 
most  refined  parts.   This  is  important  in  their  code  since  they  allow 
any  element  of  M     j  _>  2,  to  be  refined  in  creating  a  finer  space  M.. 
We  do  not,  and  the  simultaneous  displacement  iteration  considered  here 
is  perfectly  adequate  in  this  simpler  case. 

The  choice  of  smoothing  iteration  is  intimately  connected  with 
the  decomposition  of  M.  into  smooth  and  oscillatory  components.   The 
decomposition  here,  M.  =  S.  (+)  0..    is  natural  only  for  simultaneous 
displacement  iteration.   The  hard  part  of  the  analysis  here  is  the 
approximation  question  involved  in  demonstrating  the  effectiveness  of 
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step  2  of  the  multi-level  iteration  (4. 3.1)-(4. 3.3) .   For  the  decomposition 

M.  =  S.  (•+}  V.    and  the  case  considered  here  where  only  the  finest  elements 

of  M._1  are  refined  in  creating  M.,  this  approximation  question  can  be 

answered.   The  author  was  unable  to  handle  the  more  general  case  where 

any  element  of  M.  ,  may  be  refined. 
j-l 

Though  this  approximation  result  will  not  be  shown  until 
section  4.6,  we  now  state  it  since  it  is  needed  for  the  analysis  in 
this  section.   Let  {M.}.  ,  be  a  family  of  locally  refined  spaces 
satisfying  the  assumptions  of  section  4.2. 

Let  j  >_  2  and  6.  S  (0,  1)  be  arbitrary,  and  for  any 
T)  €  S .    let  T)   e  M    be  given  by 

a(n  -  n,  <f>)  =  o  ,  <(>  e  M    , 

that  is,  rj  is  the  elliptic  projection  of  n.  onto  M._  .   Then 

(4.3.4)  llln-n|||ic0el/4  !||n||,  , 

where  c  >  0  is  independent  of  j  ,  0.,  and  r\.      This  inequality,  which 

is  the  main  result  of  this  thesis,  shows  that  eigenf unctions  ty .        in 

M  whose  eigenvalues  are  relatively  small  can  be  well  approximated  in  the 

coarser  finite  element  space  M.  ,. 

J-l 

With  these  preliminaries  behind  we  are  ready  to  consider  the 
first  convergence  theorem  for  the  multi-level  iteration  (4. 3. l)-(4. 3. 3) . 
This  theorem  is  the  same  as  one  found  in  Bank  and  Dupont  for  the  case 
of  quasi-uniform  grids  and  is  similar  to  theorems  proved  by  Nicolaides, 
Hackbusch,  and  others.   Our  notation  largely  follows  Bank  and  Dupont 
although  there  are  some  differences. 
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Theorem  4.1.   There  exists  y   G  (0,  1)  and  a  positive  integer  m,  both 
independent  of  j  ^2,  such  that  if  the  approximate  solution  v  of 
(4.3.2)  satisfies 

(4.3.5)  ll!v-v!||   <y2   |||v|||    , 

then 

(<-3-6>         lllvH>n-"<J,|lliYlll«o,n-»<3)|ll  • 

Except  for  our  use  of  (4.3.4)  rather  than  a  similar  inequality  for 

quasi-uniform  grids,  the  proof  of  this  theorem  is  identical  to  the 

proof  of  the  analogous  theorem  for  quasi-uniform  grids  given  by  Bank 

and  Dupont. 

Proof.   Let  6.  £  (0,  1)  be  free  to  be  chosen  later  and  consider  the 

decomposition  of  M.: 


L   =  S.    ©  0. 

3  3     W  3 


Let  e.  ,  0  <_  £  _<  m  +  1,  be  the  convergence  error, 

e£  =   u£,n  "  u(J)  • 
since  the  second  subscript  on  u  plays  no  role  in  the  theorem.   Then 

e£  =  n£  +  k    * 
for  some  n.  €E  5  and  £  £  0.      Expanding  e_  in  eigenf unctions, 

N      .   . 

eo  =  *     al  *±   • 

i=l 

we  have  for  £  =  0,  1,...,  m, 

1-1         AN 
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by    (4.3.1).      Thus, 


gn2=    2     ai(1-^ry)2m  iii*ij)m2 

m         i=j(9)+i    x        xy}  x 

N 


Id-  e.)2m        z       a2  llli^HI2 

J  i=J(8)+l     1         1 


so 


III5JII  <  a  -  e.)m  |||5o|||  <  a  -  e.)m  |||e0|||  , 
IllnJII  i  |||n0|||  i  |||e0|||  . 

Let  £  and  ri  be  the  elliptic  projections  on  M.  .,  of  E,     and  r)  ,  respectively, 

j-1     m      m 

a(!  -  E,  (J))  =  0  ,  <j>  6  M    , 
m  j-1 

a(ii  -  n  ,  <j>)  =  0  ,  <j>  G  M    . 
m  j-1 

Since  £  is  the  projection  relative  to  the  energy  inner  product 

mi  -  gii  i  wgii  <  a  -  <>  >■  nie0in 


and  by    (4.3.4) 


n  -  njll  <  c0e1/4  Hln.HI  1  c0e1/4  |||.0| 


Since  v  =  r\  +  E,  and  e     =  n.     +  E,      , 

m  m  m 


em  -  v|||   <  c^174    |||e0|||   +  (1  -  6.)m   |||e0: 


Then  since  e    in    =  e     -  v, 
m+1  m 


llle^lll   1   ll|em-  v|||    +    [||v  -  v||| 

<cQ  \lh    |||e0|||    +   (l-8.)m    |||e0|||    +Y2    IHvlll 
by   the  hypothesis   of   the.   theorem.      Since  v  is   the  projection  of   e     on  M._-, 
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Sill  i  |||.0I 


Thus, 


(4.3-7)  Ille^JII  l(c0ej1/4+  (1  -ej)m  +  Y2)    |||e0| 

To   complete  the  proof ,    let  y  G    (0,    1)    be   such   that 

Then  choose  0.   e    (0,    1)    such   that 


c  e  <  3C 

0  j    3 


and  finally  m  such  that 


(1  -  B.)m<*  . 

Since  these  choices  of  y,  6.  and  m  are  independent  of  j,  the  proof 

is  complete.  □ 

Inequality  (4.3.5)  in  this  theorem  is  essentially  an  induction 

hypothesis.   Suppose  that  for  fixed  j   >_  2  using  iteration  (4. 3.1)- (4. 3. 3) 

with  j  =  j   causes  the  energy  norm  of  the  convergence  error  to  be  reduced 

by  the  factor  y   £  (0,  1)  according  to  this  theorem.   Also  suppose  we 

have  an  arbitrary  iteration  on  M.   .  which  reduces  the  energy  norm  of  the 

J0 
convergence  error  on  M.  _,  by  the  same  factor  y   €=  (0,  1)  at  each  iteration, 

J0 
Then  this  iteration  on  M.     can  be  used  to  approximately  solve  the 

J0 
related  problem  (4.3.2)  on  M.   -  arising  in  using  the  multi-level 

J0 
iteration  (4. 3.1)- (4. 3. 3)  on  M.  .   Beginning  with  the  trial  solution 

J0 

vQ  =  0 

.  oo 

in  M.    the  iteration  on  M.   ..  produces  a  sequence  of  iterates  {v0}„  _ 

j0-i  j0-i  *  *=0 

satisfying 

flK-v|||  </  iiivin  . 
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Thus,  two  of  these  iterations  on  M     would  be  enough  to  satisfy  the 
hypothesis  (4.3.4)  of  theorem  4.1. 

According  to  this  theorem  the  constant  y6  (0,  1)  by  which  the 
energy  norm  of  the  convergence  error  is  reduced  per  iteration  of 
(4.3.1)-(4.3.3)  is  independent  of  j  >  2.   It  follows  that  if  j  -1  >  2 
the  iteration  on  M.   .,  can  be  taken  as  the  multi-level  iteration 
(4. 3.1)-(4. 3. 3)  with  j  =  jn~l-   In  this  case  we  must  solve  in  turn 
related  problems  on  M.  ~.      Again  this  can  be  done  recursively  with 
iteration  (4.3. l)-(4. 3. 3)  if  j  -2  >  2. 

To  make  these  ideas  clear,  we  may  consider  the  following 

informal  procedure.   Calling  this  recursive  procedure  with  j  =  j 

causes  n  of  the  multi-level  iterations  (4. 3.1)-(4. 3. 3)  to  be  carried 

out  on  M.  . 
J0 


Procedure  Outer_l_on_M.  (n) 

If  (j  =  1)  {solve  the  equations  on  M.  directly;  then  return} 
else  {repeat  n  times{ 

1.  do  m  smoothing  inner  iterations  on  M.  according  to  (4.3.1). 

2.  approximately  solve  the  related  problem  (4.3.2)  on  M 

by  calling  Outer_l_on_M    (2);  then  make  the  correction  (4.3.3) 
}} 
return 


The  computational  cost  of  this  multi-level  iteration,  which 
will  be  analyzed  in  the  next  section,  depends  on  the  dimensions  of  the 
spaces  {M.}._  .   In  order  to  have  a  cost  per  iteration  proportional  to 
the  number  of  unknowns,  the  cost  of  the  smoothing  inner  iteration  on 
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M   must  dominate  the  cost  of  recursively  solving  the  related  problems  on 

j0 
the  coarser  spaces.   This  will  only  be  true  if  the  dimensions  of  the 

spaces  (M  }   ,  grow  fast  enough.   Since  these  dimensions  will  grow  quite 

slowly  for  many  families  {M ./._-.  of  locally  refined  grids  this  presents 

a  problem. 

To  ameliorate  this  difficulty,  we  consider  now  a  second 

approach  to  solving  the  related  problem  (4.3.2)  in  the  multi-level 

iteration  (4. 3. l)-(4. 3. 3) .   This  second  approach  is  based  on  a  slightly 

sharper  version  of  theorem  4.1.   First  we  require  a  simple  algebraic  lemma. 

Lemma  4.2.   Let  x,y  be  constants  in  (0,  1).   If  m  is  given  by 

[1  -  y 

m  =  *- 

I  xy 


then 


Proof.   We  have 


(1  -  x)   :  y  . 


so 


Then 


so 


..     vin+1    _ 

(1  "  X)  .   -1  =  1+  (1  -x)  +...+  (1  -x)m 
(1  -  x)  -  1 

>  (m  +  1)(1  -  x)m  , 


(1  -  x)m+1  -  1  <  -x(m  +  1)(1  -  x)m 


1  >  (1  -  x)"*1  +  x(m  +  1)(1  -  x)m  =  (1  +  mx)(l  -  x)m  , 


(1  -  x)m  <  (1  +  mx)  -1  <  (1  + Z)  1  =  y  .  D 
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Using  this  lemma,  we  can  show  the  following  theorem,  which  for 
our  purposes  improves  on  theorem  4.1.   All  notation  is  as  before. 

Theorem  4.3.   Let  r  G  (0,  1)  be  fixed  and  let  y  G  (0,  1)  and  j  >  2  be 
arbitrary.   Then  there  exists  a  positive  integer  m  independent  of  j 
such  that  if  the  approximate  solution  v  of  (4.3.2)  satisfies 

(4.3.8)  |||v  -  v|||    <   r  y    |||v|||    , 
then 

(4.3.9)  III  V!-"^  Ill  <T  lll"0-u(j)|||  . 
Moreover,  m  can  be  chosen  to  satisfy 

(4.3.10)  m  <  m  y"5 

where  m  >  0  is  a  constant  independent  of  j  and  y. 

Proof.   The  proof  is  the  same  as  that  of  theorem  4.1,  except  for  the 
choice  of  the  parameters  m,  9.,  y  at  the  end  of  the  proof.   Instead  of 
inequality  (4.3.7)  we  have 

(4.3.11)  III e^ III  1  (c^174  +  (1  -  e.)m  +  ry)  |||e0|||  , 

2 
because  the  factor  ry  in  (4.3.8)  replaces  the  factor  y   in  (4.3.5). 

Since  we  wish  to  show 

c09^/4  +  (1  -  B.)m   +  ry  <  y  , 
it  is  enough  to  choose  G .  and  m  such  that 

cne1/A  +  (i  -  e.)m  <  (i  -  r)y  . 

For  convenience  assume  c_  _>  1,  without  loss  of  generality,  and  define 


y  by 


_  (1  -  r)y 


y   =   2c0 


Then  we  may   choose  0      as 


and  m  as 
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6j        =y 


m  = 


i  -  y 


V 


Using  lemma  4.2, 

c.e}'4  +  a  -  e.)m  <   o^>l+  (l_^)l<  (1  .  r)Y  , 

establishing  the  first  part  of  the  theorem.   For  the  second  part  we  have 


m  < 


e.y 


-5 


y 


2c 


1  -  r 
2c 


0  ,5  -3 


<  1  +  (rT2-)5  Y"5 

<  (i  +  (y^)5)  y"5  , 

completing  the  proof.  D 

Using  this  theorem,  we  can  now  analyze  a  modified  version  of  the 

multi-level  iteration  given.   This  modified  iteration  will  be  the  same 

as  (4.3.1)-(4.3. 3)  except  that  m,  the  number  of  smoothing  inner  iterations, 

will  now  be  made  to  depend  on  the  level,  j.   In  applying  this  modified 

iteration  to  M.  ,  j_  >_  3,  the  related  problem  in  M,   _  will  be  solved  by 

doing  only  one  of  the  modified  outer  iterations  on  M.   _.   However,  to 

J0 
adequately  solve  this  related  problem  on  M,   -  using  only  one  outer 

J0 
iteration,  rather  than  two  as  before,  the  number  of  smoothing  inner 

iterations  must  be  increased. 
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An  informal  statement  of  this  modified  outer  iteration 
follows.   Calling  this  procedure  with  j  =  j   causes  n  multi-level 


iterations  to  be  performed  on  M.  .   Here  m(#,  •)  is  an  integer  valued 
function  to  be  given  following  the  description  of  this  procedure. 


Procedure  0uter_2_on_M.  (n,  j  ) 

If  (j  =  1)  {solve  the  equations  on  M.  directly;  then  return} 
else  {repeat  n  times{ 

1.  do  m(j ,  j  )  smoothing  inner  iterations  on  M.  according 
to  (4.3.1). 

2.  approximately  solve  the  related  problem  (4.3.2)  on  M 
by  calling  0uter_2_on_M.    (1,  j  ),  then  make  the 
correction  (4.3.3). 

}} 
return 


It  remains  only  to  choose  the  function  m(* ,  •)  required  by 
this  procedure.   Let  r  €E  (0,  1)  be  arbitrary  and  let  m  >  0  be  the 
corresponding  constant  from  inequality  (4.3.10)  of  theorem  4.3.   Also,  let 
Yn  e  (0,  1)  be  arbitrary.   Then  we  may  set 


(4.3.12)  m(j,  jQ)  = 


-5   5«-V 


mo  yo    r 


Since  the  argument  needed  to  show  that  the  multi-level  iteration  in 

0uter_2  converges  rapidly  is  somewhat  more  complex  than  that  for 

r     1°° 
Outer  1,  we  state  this  as  a  theorem.   Let  iu^  j      .,  be  successive 
—  0,n  n=l 

iterates  produced  by  the  multi-level  iteration  in  0uter_2.   Then  we  have 
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Theorem  4.4.   Let  J   >  2  be  arbitrary.   Then  with  m(*,  •)  given  by 

(4.3.12)  the  iterates  {u.   }   ,  satisfy 

0,n  n=l 

(jn)  (jn) 

I^W  °  "I  £T0  UI-o.n  "  "     I"'  - 

for  all  n  >_  1. 

Proof.   In  order  for  each  outer  iteration  on  M.  ,  i.  >  2,  to  reduce  the 

energy  norm  of  the  convergence  error  by  the  factor  y    ,    according  to 

theorem  4.3,  it  suffices  to  do  m  of  the  inner  iterations  (4.3.1)  on 

M.   for  some 

m  <  mQ  Yq5  , 

and  then  approximately  solve  the  related  problem  (4.3.2)  with  relative 
accuracy  ry  .   Examination  of  (4.3.11)  shows  that  doing  extra  inner 
iterations  never  hurts  so  the  choice 

-5 


m  =  m(J0>  3Q)   = 


Lm0Y0 


is  fine.   Thus,  for  the  case  where  j_  =  2  and  the  related  problem  on 

M.   .  is  solved  directly  the  proof  is  done.   Otherwise,  we  solve  the 

30 
related  problem  on  M.    with  one  outer  iteration.   Choosing  Y  =  rY„ 

30 
in  theorem  4.3  it  suffices  to  do  m  inner  iterations  on  M.   n  for  some 

Jo"1 

m  <   mQ(rY0) 

and  then  solve  the  related  problem  on  M.   0  with  relative  accuracy 
r  YQ.   Doing 

/•   -,     n    I     -5-5 
m  =  m(j0-l,  j0)  =  LmQ  YQ  r 

inner  iterations  is  certainly  sufficient,  so  the  proof  is  done  for 

jn  =  3.   Continuing  the  induction  in  the  obvious  way  completes  the  proof.  D 
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The  computational  cost  of  the  multi-level  iterations  described 
in  this  section  and  construction  of  an  attractive  overall  algorithm 
for  obtaining  solutions  to  the  elliptic  problem  (2.1.3)  will  be 
considered  in  the  following  section. 

4.4.   Complexity 

In  this  section  an  estimate  of  the  computational  cost  of  the 

multi-level  iterations  described  in  the  last  section  is  given,  and  these 

iterations  are  then  used  as  building  blocks  for  computationally  attractive 

algorithms  for  solving  finite  element  equations  on  locally  refined  grids. 

Under  certain  restrictions  these  algorithms  will  be  shown  to  be  of 

optimal  order,  producing  a  good  solution  to  the  finite  element  equations 

for  M   k  2l  1,  in  0(N,  )  operations,  where  N  is  the  dimension  of  the 

finite  element  space  M,  .   This  complexity  analysis  is  somewhat  harder 

than  the  corresponding  analysis  in  chapter  three,  since  the  complexity 

of  the  multi-level  iterations  in  the  last  section  depends  on  the 

dimensions  of  the  finite  element  spaces  {M.}.  ..  used.   This  complication 

3    j=l 

was  absent  from  the  last  chapter,  since  for  the  family  {M.}.  ,  of 

J  J=l 

quasi-uniform  spaces  there,  the  dimensions  of  the  spaces  were  essentially 

determined  by  the  dimension  of  the  coarsest  space  M  and  the  mesh  ratio  p . 

A  second  minor  difficulty  also  arises  in  proving  the  optimal 

order  convergence  of  the  algorithms  here.   A  multi-level  algorithm  is 

considered  to  be  of  optimal  order  if,  when  applied  to  any  finite  element 

space  M,  in  {M.}._,,  it  produces  a  good  solution  u    in  $(N,)  operations, 
k      J  j— 1  k 

To  define  what  "a  good  solution"  means,  an  error  bound  applicable  to 
locally  refined  grids  is  needed.   We  now  proceed  to  derive  such  a  bound. 
It  is  important  here  that  the  error  bound  derived  not  be  overly 
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pessimistic.   If  it  were,  the  results  here  would  be  correspondingly 

weakened.   To  have  shown  a  multi-level  method  produces  a  solution 

containing  a  convergence  error  smaller  than  the  expected  truncation 

error  means  very  little  if  one  expects  far  more  truncation  error  than 

actually  results. 

Only  a  simple  error  estimate  is  required  here.   Such  an 

estimate  can  be  derived  directly  from  the  Bramble-Hilbert  lemma.   First 

we  need  the  notion  of  a  Sobolev  seminorm.   For  I    a  non-negative  integer, 

the  seminorm  |  •  |  .is  the  same  as  the  Sobolev  norm  ||  •  ||  .  except  that  it 

IT  if 

contains  only  derivatives  of  order  £,  not  those  of  lower  order.   That 

is,  for  u  £  H  (fi) , 

|u|      =  (  I      l^ull2/'2  . 

h*(B)        |o|-£        l 

Just  as  for  the  Sobolev  norm,  these  seminorms  can  be  defined  for 

non- integral  I   in  a  more  complicated  way. 

Now  let  M  be  a  finite  element  space  of  C  Lagrangian  elements, 

based  on  a  triangulation  T  of  Q   and  satisfying  conditions  (4. 2. 1)- (4. 2. 5) 

Suppose  the  element  trial  functions  are  polynomials  of  degree  k  -  1. 

Then  M  is  said  to  be  a  finite  element  space  of  degree  k.   For 

u  G  H(J2),  1  <  I  <   k,  the  nodal  interpolant  u  G  M  of  u  is  well  defined 
E  1 

and  we  have 


Lemma  4.5.   For  any  triangle  T  £  T  and  for  s  =  0,  1 

U  "  UTI  s    -  as  dT    lUl  I 

1  HS(T)     S  L  H^CT) 

where  d   is  the  diameter  of  T. 
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The  variable  i   here  need  not  be  an  integer.   This  lemma  is  a  special 
case  of  the  Bramble-Hilbert  lemma,  which  may  be  found  in,  for  example, 
Strong  and  Fix  or  Oden  and  Reddy. 

r    i  °° 

Applying  this  lemma  to  the  spaces  in  a  family  {M.}.  ,  of 

3   J=1 

locally  refined  finite  element  spaces  satisfying  the  assumptions  of 

section  4.2  yields  the  required  error  bound.   Assume  that  the 

1+a 
differential  equation  (2.1.3)  is  H    regular  for  a^  (0,  1),  that 

1+a 
is,  that  its  solutions  lie  in  H   .   Then  the  following  error  estimate 

holds  on  M. ,  i  >  1. 
J    - 


Theorem  4.6.   Let  u    be  the  finite  element  solution  of  (2.1.3)  in 

M . .   Then 
J 

(4.4.1)  |||u-u^|||    <cq(E        df    |u|21+a        )1/2 

U  T^T.  L  H1  a(T) 

where  c   >  0  is  independent  of  u,  j,  and  the  family  of  locally  refined 


00 

spaces  {M . } . 

i  i=l 


Proof.      Let   u  be   the   interpolant   of   u   in  M..      Then 

III-   -   u(j)|||    =    min     |||u-   v|||    <    |||u-u<j) 
vSM. 

Since  the  energy  and  H  norms  are  equivalent, 

(4.4.2)  l||u-u(j)|||  <  c  Mu-u^M    . 

H 


Applying  the  above  lemma  on  each  triangle  in  T . , 

M      (j)M2      .2       v   ,2  (1+a)  ,  .2 
I  u  "  V    2    -  an   Z   dT       u  l+« 

1   L2(fi)     °^T.   T         H1+a(T) 


Thus, 
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(j)l2     ^  J-       v   ^2a  I  I2 

U   "   UT     I   1         <   °1      ^     d        U     1_L~ 


(j)ii2    .  2     _  .  ,2a  ,  .2 
u  "  u tJ  II  !    <  c    E   d   |u| 

HX(fi)      T^T.        H±W(T) 


2        2   2 
for  c  =  max{an,  a  } ,  assuming  d  <  1  for  all  triangles  T  £  T. 


Now 
J 

combining  this  with  (A. 4. 2)  the  result  follows.  D 


Since  this  bound  contains  d  ,  in  effect  the  local  mesh  parameter, 
it  is  much  tighter  for  locally  refined  grids  than  bounds  in  terms  of  the 
maximum  element  diameter.   Moreover,  the  exponent  of  d  here  is  optimal. 
For  sufficiently  smooth  u  the  reverse  inequality,  bounding  the  error 
from  below,  may  be  expected  to  hold  although  this  is  difficult  to  show. 
Establishing  that  the  multi-level  methods  of  this  section  produce 
solutions  whose  convergence  error  is  comparable  to  the  truncation  error 
expected  from  this  error  bound  may  be  regarded  as  more  than  adequate  for 
most  purposes.   In  particular,  for  problems  with  singular  solutions 
the  improved  convergence  of  finite  element  spaces  using  the  proper 
locally  refined  grids  may  be  readily  shown  using  this  bound.   The 
multi-level  algorithms  here  will  inherit  this  improved  convergence  to 
singular  solutions  on  properly  refined  grids. 

With  this  error  bound  established  we  are  ready  to  turn  to  the 
design  of  multi-level  algorithms  for  solving  the  elliptic  problem  (2.1.3) 

00 

on  locally  refined  grids.   Let  IM.}._1  be  a  family  of  locally  refined 

finite  element  spaces  satisfying  the  assumptions  of  section  4.2.   The 

~(k) 
algorithms  we  wish  to  consider  will  produce  a  good  solution  u    in  any 

number  M  of  this  family,  k  >_  1.   The  first  step  in  these  algorithms  is 
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to  solve  the  finite  element  equations  for  M,  directly  to  obtain  an 
approximate  solution  u    to  the  problem  (2.1.3).   After  this,  solutions 
u    £  M,,  2  <  j  <  k,  will  be  produced,  each  a  good  approximate  solution 
of  (2.1.3).   Each  solution  u    ,  2  <_  j  <_  k,  is  obtained  by  beginning 
with  the  previous  solution  u      G  M    ,  considering  it  instead  as  a 
function  in  M . ,  and  then  carrying  out  a  sequence  of  multi-level  iterations 
to  obtain  the  improved  solution  u    G  M..   This  multi-level  iteration 
may  be  taken  as  either  of  the  algorithms  of  the  last  section. 

Two  questions  arise  immediately:   what  conditions  must  the 
above  process  satisfy  in  order  to  obtain  an  accurate  solution  u    , 
and  what  determines  the  computational  complexity  of  this  process?  With 

the  error  bound  (4.4.1)  already  proven,  the  first  question  is  easier. 

~  (k) 
We  wish  to  obtain  a  solution  u    which  contains  only  convergence  error 

comparable  to  the  expected  truncation  error.   That  is,  we  want  the 

convergence  error  to  satisfy 

III  00     (k)|ll  •    t    v  a2ql    I  I2      \1/2 
WW         -  «   III  1  c  (  E   d    |u|       ) 

u  T^T,  l  H  ™(T) 

k 

~(k) 
In  order  to  have  an  algorithm  whose  solution  u    satisfies  this 

inequality,  it  is  convenient  to  ask  that  the  intermediate  solutions 

u   ,  1  j<  j  £  k  -  1,  also  satisfy  this  type  of  inequaltiy.   That  is, 

for  1  <   j  <   k, 

(4.4.3)      |||5(J)  -u^HI  <c0(  I   df  |u|21+a   )1/2  . 

U  T^T.  L  H  ^(T) 

3 

A  simple  condition  insuring  this  is  given  in  the  following  theorem. 
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Theorem  4.7.   Suppose  that  u*   satisfies  (A. 4. 3),  and  for  2   <   j  <   k 

(4.4.4)     ma"'  -»«>|||  i  a  +  21^)-1  iBstt-"  -  .«>«  . 

Then  (4.4.3)  holds  for  1  <   j  <  k. 

Proof.   By  assumption 

IH~(1)    (Din  .        f    -   ,2a  |  ,2      ,1/2 
llu    "  u    III  1  c  (  Z   d    u        ) 

0  reTl  T    H1+a(T) 

Now  inductively  suppose 

|||~(j)    (j)ii|  <        f    v       ,2a  |  ,2      .1/2 
III11    "  u    III  1  cn(  Z   dT   u  1-Kv   ' 

0  TET   T     H1+a(T) 

for  some  j  ,  1  <_  j  <^  k.   Then 

|||5(j)  -.|||  <  |||5(3)  -u(3)|||  +  |||u(j)  -.|||  ; 
so 

|||5«>  -  ulH  <  2c0(  E   df  |u|2^   )1/2 
0  rer.  T    h14<*(t) 

using  the  error  bound  (4.4.1).   Then,  since 

III    (J+1)  III    *        t     v  A2a  I     I2  ^/2 

llu  -  ulll   1  cn(     E  dT       u     14a        ) 

°   TGT..-       T  H14a(T) 

3+1 

we  have 

(4.4.5)     |||u^.U^+1>|||<c0(  E    df  |u|21+a   )1/2 

0  r^T.+1  T    H1+a(T) 

+  2c  (  I      d    u        ) 
0  Tex.  T    h14<*(t) 

Now  let  T  6  T.  and  let  S  be  the  set  of  subtriangles  of  T  in  T.,-,. 
That  is  S   consists  either  of  T  or  of  the  two  or  four  triangles  formed 
by  refining  T.   Then  we  have  the  following  relations. 
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dT<2dT,  ,  ve   sT 

II2  v         I  I2 

,u|  1+a    "   Z   ,u|  1+a 

H  a(T)    T'eST    H  ™(T')  . 

Applying  these  to  (4.4.5)  yields 

|||~(j)    (J+1)||l  <        (     y  ,2a  i  ,2      .1/2 

lllu    -  u     HI  <  c  (   I    d    |u|       ) 

T^T  H    (T) 

+  2c  (  Z     E    (2d  )2a  |u|2       )1/2 
U  T6T.  T'eST  H  ^(T') 

=  cn(l  +  2(2)a)  (   Z     df  |u|         )1/2  . 
U  T^Tj+1   L     H1+a(T') 

Finally,  making  use  of  the  hypothesis  (4.4.4)  the  induction  is  complete.   D 


The  dependence  of  this  result  on  the  regularity  of  the  problem 
will  affect  the  complexity  of  the  multi-level  algorithm  here  in  an 
interesting  way.   Consider  the  construction  of  such  algorithms  based  on 
this  theorem.   As  described  above,  after  obtaining  u  "   by  solving  the 
equations  for  M   directly,  approximate  solutions  u    eM.,  2j£j<^k, 
are  obtained  by  iteratively  improving  the  previous  solution  u      €=  ^-_i' 
The  amount  of  iterative  improvement  required  is  given  by  (4.4.5).   If  the 
problem  is  H    regular,  we  must  reduce  the  convergence  error  by  the 
factor  (1+2   )     in  going  from  u      to  u    .   The  interesting  point 
here  is  that  more  error  reduction  is  required  for  smoother  problems. 
Since  the  regularity  did  not  enter  the  convergence  theorems  of  the  last 
section,  obtaining  a  given  amount  of  reduction  of  the  convergence  error 
is  apparently  no  easier  for  problems  with  smooth  solutions  than  for 
problems  with  weak  regularity.   Thus,  the  computational  cost  per  element 
of  the  algorithms  here  will  be  greater  for  problems  with  smooth  solutions, 
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Of  course,  there  is  no  real  dilemma  here.   The  algorithms 

here  produce  a  solution  satisfying  (4.4.5).   Since  the  finite  element 

(k) 
solution  u    satisfies  (4.4.1),  the  approximate  solution  produced  by 

the  multi-level  algorithm  must  satisfy 

|||u<k>  -  u||!  <2c.(  I   df  \u\2  )l>2    • 

a  much  stronger  inequality  for  large  a.   Thus,  even  though  one  must  work 

harder  per  element  for  smooth  problems,  it  will  still  be  cheaper  to 

produce  a  good  answer. 

Now  we  are  ready  to  consider  the  design  of  multi-level 

algorithms  based  on  theorem  4.7  and  the  multi- level  iterations  of 

section  4.3.   The  multi-level  iteration  used  may  be  that  described  in 

either  procedure  Outer_l  or  0uter_2.   In  either  case  the  dimensions  of 

the  finite  element  space  {M .  s .    -,  will  be  of  central  importance.   If  N. 

is  the  dimension  of  M.,  j  _>  1,  optimal  complexity  bounds  can  be  shown 

under  the  assumption  that  the  dimensions  {N.}.  .  grow  at  least 

3  1=1 

geometrically.   That  is, 

(4.4.6)  N...  '  <5N.  ,  j  >  1  , 

for  fixed  6  >  1.   Actually,  this  can  be  weakened  slightly.   It  is  enough 

oo 

that  we  have  a  sequence  of  numbers  {x,}._   growing  at  least  geometrically, 

*j+1  1  6x.  ,  j  >  1  , 

for  6  >  1,  and  that  N.  and  x.  are  comparable, 

J      J 

—  x.  <  N .  <  c  x. 
c  j  —  3  -    J 

for  fixed  c  >  0.   For  simplicity,  assume  that  (4.4.6)  holds.   Using  either 
of  the  multi-level  iterations  of  section  4.3,  the  storage  required  to 
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approximately  solve  the  finite  element  equations  for  M  ,  k  _>  1,  can  be 
estimated  as  follows.   Let  c  be  the  storage  required  for  direct 
solutions  of  the  equations  for  M  .   Then  the  total  storage  required 
is  proportional  to 


k  k 

I  N.  +  c,  <  N.   Z   6J   +  c, 

3=2  3    ^  kj  =  2         1 

<  \  t^t  +  Cl  • 

Since  c   is  fixed,  for  any  6  >  1  the  storage  required  is  proportional 

to  N,  . 
k 

Next  consider  the  use  of  the  multi-level  iteration  described 
by  procedure  Outer_l.   According  to  theorem  4.1,  when  this  multi-level 
iteration  is  applied  to  M.,  j  >  1,  the  convergence  error  is  reduced  by 
a  factor  Y  at  every  iteration,  where  y  <   1  is  independent  of  j. 
Actually,  the  proof  of  theorem  4.1  shows  Y  can  be  chosen  arbitrarily, 
though  this  is  unimportant,  since  for  any  y  <   1  ~we   can  find  V  ^>  1 
such  that 

YV  <  (1  +  21+V1  . 

Let  y  <  1  and  an  integer  v  satisfying  this  inequality  be  fixed.   Then 

~(k) 
according  to  theorem  4.7,  a  solution  u    in  M  ,  satisfying  (4.4.5)  can 

K. 

be  produced  by  carrying  out  V  of  the  multi-level  iterations  described 
by  Outer_l  on  each  space  M.  ,  2  _<  j  <^  k.   For  each  M.  ,  2  j<  j  j<  k,  the 
approximate  solution  produced  in  the  next  coarser  space  M    is  used 
as  a  starting  value  for  this  iteration.   To  estimate  the  complexity  of 
this  procedure,  note  that  in  applying  Outer_l  to  M.,  j  >_  2,  0(N.) 
operations  are  required  plus  the  cost  of  recursively  solving  the  related 
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problem  in  M   , .   This  related  problem  is  solved  in  Outer_l  applied  to 
M.,  by  recursively  applying  Outer_l  twice  to  the  equations  for  M   , ,  if 
j  -  1  >  1.   Continuing  the  recursion  one  sees  that  the  total  cost  in 
doing  V  multi-level  iterations  of  this  type  on  the  equations  for  M.  is 
proportional  to 

v(  Z   2j_i  N  +  2j_1  cj  <  VN.(  Z   (|)j_i  +  (|)j_1  cj  , 

i=2       1         l     ~       J  i=2   6        6      2 

where  c„  is  the  cost  of  solving  the  linear  system  for  M  .   Since  c„  is 

a  fixed  constant,  the  computation  time  T(M.)  to  do  v  of  these  multi- 

J 

level  iterations  on  M.  must  satisfy 

f 

0(N.)  ,  6  >  2 

T(M.)  <  J   0(N.  log  (N.))  ,   6=2 

log  2 
L0(N.  log  6)  ,     1  <  6  <  2  . 

Then  following  the  procedure  described  above  where  the  approximate  solution 

~(k) 

u    in  M  is  obtained  by  successively  obtaining  good  solutions  in 

M.,  1  <^  j  <_  k,  the  total  computation  time  T  to  obtain  u    must  satisfy 

(\) 

T  =  Z  T(M.)  <  I   (N,  log  (N,  ))  ,6  =  2 
j=l    J      \       k 

log  2 

[(Nk  l0g  6)  ,     1  <  6  <  2  . 

For  many  problems  a  bound  like  this  is  quite  adequate.   For 

example,  in  doing  local  refinement  along  a  plane  singularity  or  boundery 

in  a  three  dimensional  problem,  the  dimensions  {N.}.  ,  would  at  least 

J  J=l 

quadruple  at  each  level.   Since  the  two  dimensional  results  here 
readily  extend  to  three  dimensional  problems,  this  example  is  important. 
In  other  problems,  for  example,  refining  along  a  line  sigularity  in  two 
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or  three  dimensions,  the  dimensions  typically  only  double  at  each  level, 
yielding  the  N  log  (N)  bound.   Still  other  problems  arise  for  which  the 
growth  of  these  dimensions  is  even  slower.   In  refining  around  point 
singularities,  these  dimensions  will  usually  not  grow  geometrically, 
and  the  results  of  this  section  are  of  little  use.   This  does  not 
necessarily  mean  that  multi-level  methods  should  be  avoided  in  these 
cases,  only  that  theoretical  results  showing  their  effectiveness  have 
not  yet  been  given. 

When  the  dimensions  of  the  spaces  M.  grow  geometrically 
according  to  (A. 4. 7),  but  the  rate  of  growth,  5,  is  not  greater  than 
2,  the  complexity  results  above  can  be  improved  by  replacing  the 
multi-level  iteration  given  by  Outer_l  with  that  given  by  0uter_2. 
In  equation  (4.3.12),  governing  the  number  of  inner  iterations  used  on 
each  level  in  0uter_2,  we  are  free  to  choose  both  yn  and  r.   For 
convenience,  yn   may  be  chosen  as 


With  this  choice  in  applying  0uter_2  to  the  equation  for  M^  , 

0 


J0(j  ) 

2  <_   j   <_  k,  only  one  iteration  of  0uter_2  is  needed  to  get  u  0 

satisfying  (4.4.4)  from  u      .   The  other  free  parameter,  r,  in 


(4.3.12)  will  be  chosen  as 


r  =  6   u  . 
Then  for  any  6  >  1,  r  6  (0,  1),  so  theorem  4.4  shows  that  applying  the 
multi-level  iteration  described  by  Outer  2  to  M.  ,  2  <  i.  <  k,  reduces 

V        °" 

the  convergence  error  by  Yn  Per  iteration  as  required.   From  equation 
(4.3.12),  governing  the  number  of  inner  iterations  on  each  level,  we  can 
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~(V    ~(V1} 

estimate  the  computation  time  to  produce  u     from  u      ,  j   >  2. 
This  time,  T(M.  )  is  proportional  to 

j0  -5(jn-j)  j0    -5  j0_j 

E   r         N  +  c  <  N    E   (f  )     +  c. 
j=2  J    Z    J0  j-2  5  2 

j0    -I  .   . 

<  N.   E   (6  6)j0  j  +  c 
J0  j-2  X 


<  N.   (1-6   6)  1  +  c9  . 

Thus,  the  amount  of  computation  required  to  produce  u     from  u 

j   >  2,  is  0(N.  ).   Then,  since 
J0 


k  -1  -1 

E   N.  <   N   (1  -  6  X) 

j0=2  ^0 

for  any  6  >  1,  the  total  computation  time  is  0(N  ). 

This  is  an  optimal  bound,  showing  0(N)  complexity  under  the 
weaker  assumption  that  the  dimensions  IN.}.,  grow  geometrically,  with 
any  growth  rate  6  >  1.   Unfortunately,  the  dependence  of  this  0(N)  bound 
on  6  is  quite  severe,  though  this  could  be  improved  by  choosing  the 
various  constants  differently,  and  by  using  a  more  effective  inner 
iteration.   The  hidden  constant  in  the  0(N)  bound  just  given  contains 
not  only  the  factors  (1-6    )    and  (1-6   )    seen  above,  but 
also  the  constant  m  of  (4.3.12),  which  depends  indirectly  on  6.   By 
the  proof  of  theorem  4.3  it  can  be  seen  that  mn  goes  as 
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(1  -  r)  5  =  (1  -  6   6)  5 

Thus,  including  the  dependence  on  6,  the  complexity  bound  is 

1 
0((1  -  6  V6-  (1  -  6"1)"1  N)  . 

_1  /£  _C.  _-|   _  1  C. 

While  the  factor  (1-6    )    (1  -  6   )    already  exceeds  10   for 
6=2,  in  practice  one  would  expect  the  dependence  on  6  to  be  far 
less  severe  than  this,  although  the  appropriate  numerical  experiments 
have  not  yet  been  carried  out. 

This  completes  the  analysis  of  the  algebraic  aspects  of  the 
multi-level  algorithms  of  this  chapter.   It  remains  to  prove  the 
approximation  result  (4.3.4),  on  which  the  complexity  results  just 
given  are  based. 

4.5.   Interpolation 

This  section  begins  the  analysis  of  the  approximation  properties 
of  locally  refined  finite  element  spaces  needed  to  justify  the  convergence 
and  complexity  results  of  the  last  two  sections.   For  convenience  this 
analysis  is  split  into  two  parts.   The  first  of  these,  given  in  this 
section,  deals  with  properties  of  the  interpolation  mapping  from  M.  to 

00 

M.  ,,  i  >  2,  for  M.,  M.  ,  members  of  the  family  {M.}.  ,  of  locally 
J-l    -         3    3-1  3  3=1 

refined  spaces.   The  second  part  of  this  analysis,  given  in  the 
following  section  will  use  the  results  of  this  section  to  obtain  the 
required  approximation  theorem. 

While  the  analysis  in  the  next  section  is  completely  general, 
the  analysis  here  is  specific  to  the  C   triangular  elements  and  grids 
based  on  aligned  triangulations  being  considered.   However,  there  is  no 
reason  to  believe  that  similar  results  could  not  be  shown  for  more  general 
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grids  and  for  most  common  elements  including  triangular  elements  with 
higher  continuity,  square  or  quadrilateral  elements,  and  elements  with 
curved  edges.   The  results  here  were  not  proven  in  greater  generality 
primarily  because  it  did  not  seem  to  aid  understanding  of  the  multi-level 
algorithms  here  and  because  the  notations  and  arguments  required  would 
have  been  burdensome  to  both  author  and  reader.   Proof  of  results 
analogous  to  those  in  this  section  for  most  reasonable  finite  element 
spaces  is  not  difficult. 

The  work  in  this  section  is  closely  related  to  the  work  Bank 
and  Dupont  did  in  analyzing  their  two-level  scheme  discussed  briefly  in 
section  4.1.   Much  of  the  notation  here  and  one  of  the  lemmas,  lemma  4.10, 
are  taken  directly  from  their  paper.   The  difference  here  is  mainly  in 
emphasis.   While  lemma  4.10  plays  a  fairly  small  role  here,  it  is 
central  to  the  convergence  proof  for  the  two-level  scheme  in  their 
paper.   Lemma  4.11,  which  is  not  from  their  paper,  is  really  much  more 
central  here. 

The  reason  for  considering  first  the  interpolation  mapping  is 
the  usual  one:   while  projection  with  respect  to  the  energy  inner  product 
is  global,  interpolation  gives  a  local  projection  more  easily  analyzed. 
Yet  the  two  mappings  are  sufficiently  alike  that  properties  of  one  can  be 
derived  from  the  other.   This  is  a  common  point  of  view,  but  the  way  in 
which  properties  of  the  elliptic  projection  are  derived  here  from  the 
interpolation  projection  is  unusual.   Instead  of  bounding  the  approxi- 
mation error  in  the  interpolate  and  using  this  as  a  bound  on  the 
approximation  error  in  the  elliptic  projection,  we  begin  by  establishing 
properties  of  the  null  space  of  the  interpolation  mapping.   From  these 
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properties,  in  the  next  section  we  will  derive  properties  of  the  null 

space  of  the  energy  projection.   Finally  these  lead  to  the  required 

approximation  result. 

This  proof  may  seem  circuitous,  but  the  author  was  unable  to 

establish  an  approximation  result  such  as  (4.3.4)  in  any  other  way. 

The  problem  encountered  with  the  usual  finite  element  approach  is  not 

hard  to  see.   Since  elliptic  projection  is  global,  an  error  estimate, 

such  as  theorem  4.6,  for  a  locally  refined  space  M,  must  involve  the 

maximum  element  diameter  h    .   However,  the  eigenvalue  bounds  involve 

max 

h  .  .   Thus,  one  will  be  able  to  prove  inequalities  such  as  (4.3.4)  but 

mm  . 

h 

the  constant  c  occuring  will  contain  some  power  of  the  factor  — . 

min 

This  factor  becomes  unbounded  for  the  locally  refined  grids  being 

considered,  so  such  a  result  would  be  too  weak  to  prove  0(N) 

complexity  bounds  like  those  in  the  last  section. 

Since  the  finite  element  spaces  considered  in  this  chapter 

are  based  on  C  Lagrangian  elements,  the  interpolation  mapping  from  M. 

to  M.  , ,  j  >  2,  with  M.,  M.  ,  members  of  {M.}.  ,,  is  automatically  well 
J-l    -  J    3-1  3  3=1 

defined.   For  any  u  e  M . ,  u  is  evaluated  at  the  points  corresponding  to 

the  nodal  parameters  of  M.  ...   Then  there  is  a  unique  uT  €=  M.  .  attaining 

3-1  I    3-1 

the  same  values  at  these  points.   One  of  the  problems  encountered  when 

using  elements  with  higher  continuity  is  that  this  interpolation  may  not 

be  well  defined.   There  are,  for  example,  C   elements  using  first 

derivatives  as  nodal  parameters  and  almost  all  C  elements  use  second 

derivatives.   In  these  cases  care  must  be  taken  that  the  nodes  of  M.  .. 

J-l 

lie  at  points  where  the  corresponding  values  of  functions  in  M.  are 
well  defined.   In  those  rare  cases  where  nodal  interpolation  is  not  well 
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defined,  more  complex  interpolation  projections  based  on  various  local 

averaging  techniques  must  be  used. 

Now  let  M  ,  j  _>  2,  denote  the  null  space  of  the  interpolation 

from  M,  to  M,  ..   M.  is  the  subspace  of  M,  consisting  of  functions 

vanishing  at  the  nodes  of  M.  ..   Much  of  this  section  and  part  of  the 

j-l 

next  are  devoted  to  establishing  properties  of  these  null  spaces.   For 
now  we  note  that  M.  yields  a  direct  sum  decomposition  of  M.: 

J    J  w  j-l 

This  is  easy  to  show,  for  example,  by  dimensionality  arguments. 

Besides  the  direct  method  of  defining  the  interpolate 

uT  G  M .    of  a  function  u  G  M.,  there  is  also  a  more  abstract  way  which 
I    J-l  J 

we  now  consider. 

The  interpolation  mapping  from  M.  to  M,_   considered  here 

is  completely  local  in  the  sense  that  the  value  of  the  interpolate 

u,.  £  M.  ,  on  a  triangle  T  6  T.  ,  depends  only  on  the  value  of  u  £  M. 
I    J-l  J-l  J 

on  T.   For  this  reason  the  interpolation  mapping  from  M.  to  M.  , 

J     J-l 

decomposes  into  separate  mappings  from  M. |   to  Ml   for  each 
T  G  T-_-i*   Rather  than  considering  separately  each  of  the  function 
spaces  M.|   and  M    |T  as  T  ranges  over  M.   we  can  consider  only  a 
few  canonical  function  spaces  on  T_  and  the  interpolation  mappings 
between  these  canonical  spaces.   These  canonical  spaces  on  T^  can  then 
be  transformed  into  the  spaces  M. |   and  M    |   through  the  affine 
transformation  F  . 

As  in  section  4.2  let  L   be  the  element  trial  space  on  the 
reference  triangle  T  .   We  must  consider  three  different  refinements 

is. 

of  T_:   regular  refinement,  green  refinement,  and  the  null  refinement, 
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in  which  T_  is  not  altered.   There  are  three  possible  green  refinements 
R 

but  we  need  only  one,  say  the  one  with  the  new  edge  leading  to  the  top 
vertex. 


Regular  Green  Null 

Figure  16.  Refinements  of  T 

R 

Since  each  of  these  refinements  of  TD  is  a  triangulation, 

finite  element  spaces  on  each  of  them  can  be  generated  in  the  usual  way 

from  the  reference  element  trial  space  L. 

Let  K     be  the  finite  element  space  on  T   generated  in  this 
R  R 

way  from  the  regular  refinement  of  T  and  let  Kr   be  the  finite  element 

R  Vj 

space  generated  by  this  green  refinement.   Null  refinement  just 

reproduces  L.      All  three  of  these  spaces  generated  by  refinement,  K  , 

R 

Kc ,  and  L   contain  L   as  a  subspace.   Moreover,  there  are  interpolation 
mappings  I_ ,  I„ ,  and  I.  from  K   ,  Kn,    and  L,    respectively  to  L.      The 

R    U         L         R    (jt 

last  is  of  course  just  the  identity  on  L.      Just  as  for  the  interpolation 

mapping  from  M.  to  M.   ,  we  denote  the  kernels  of  I.,  and  I„  by  K.     and 
j      j-1  R      b     R 

Kr ,  respectively.   The  kernel  of  I.  is  the  trivial  function  space  {0}. 

Now  let  T  G  T.  ,  be  arbitrary.   In  going  from  T.  n  to  T,  the 
J-1  .1-1     J 

triangle  T  is  regular,  green,  or  null  refined.   Consider  the  same  type 

of  refinement  of  T  .   Then  there  is  an  affine  transformation 

R 

F  :  T  ->  T  carrying  the  subtriangles  of  T0  generated  by  this  refinement 
1    R  R 

to  the  subtriangles  created  in  refining  T.   When  T  is  regular  or  null 
refined,  any  of  the  six  affine  transformations  have  this  property.   When 
T  is  green  refined  two  of  the  six  affine  transformations  carry  the  green 
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subtriangles  of  T_  to  those  of  T.   The  other  four  affine  transformations 

K 

correspond  to  the  other  possible  green  refinements  of  T. 

To  see  the  importance  of  all  this,  let  I  be  the  interpolation 

mapping  from  ML  to  M   . |   and  let  F'  be  the  dual  of  F  ,  carrying 

functions  on  T  to  functions  on  T   , 

R 


by 


F^,:  H^T)  ■*  H1^) 


F.J,(u)  =  u  o  FT  . 


For  each  of  the  three  types  of  refinement  of  T  we  get  a  commutative 
diagram.   For  regular  refinement  we  have 


F' 


L  '  T    M.  ,1 

J-1'T 

tXR  tXT 


T 
R  "j'T 


K  «— ^ —   M 


For  green  refinement  we  have 


F' 
L  h i M 


f< 


1-1 'T 
t1! 


F' 
KG  —^       WllT 


And  finally  for  null  refinement  we  have 
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f; 

L  •* —  M.  ,1 

j-l'T 

■  FT 

L  ^-^-M.|T 

That  these  diagrams  are  in  fact  commutative  is  a  consequence 

of  the  fact  that  F   carries  not  only  T   and  its  subtriangles  into  T  and 

its  subtriangles,  but  carries  nodes  of  T   and  its  subtriangles  into 

R 

corresponding  nodes  of  T  and  its  subtriangles.   We  omit  the  straight- 
forward but  messy  verification. 

The  lemmas  which  follow  just  elucidate  the  properties  of 
these  three  commutative  diagrams  and  the  various  function  spaces  involved, 

The  first  result  here  shows  that  functions  in  L   cannot  well  approximate 

•\     ^\ 
nonzero  functions  in  K.     or  K  ,  the  kernels  of  I  and  I  .   The  kernel, 

R     G  R      G 

{0},  of  I.  contains  no  nonzero  functions  so  the  question  does  not  arise 

/\      ^ 
there.   The  functions  in  K_  and  Kn   vanish  at  the  nodes  of  T-  while 

R       G  R 

nonzero  functions  in  L   cannot  vanish  at  all  nodes  of  T_.   The  results 

R 

follow  from  this  and  the  finite  dimensionality  of  all  spaces  involved. 
For  convenience  define  the  interpolation 

nL:    C°(TR)  ->  L    . 

The  mappings  I_ ,  In ,  and  I.  described  above  are  just  restrictions  of 
R    G        L 

7T.  to  the  subspaces  K  ,  K    ,    and  L   of  C  (T  ) ,  respectively. 

L  R    G  R 


Lemma  4.8.   If  w  £  K  or  w  6  K  ,  for  any  v  £  L 

R 
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(4.5.1a)  l|w-v||  >c||w|| 

LZ(TR)  l/(TR) 

(A. 5.1b)  |w  -   v|     .  >_  c    |  w  |    ..  , 

El(TR)  '  tiLCTR) 

for  c   >   0   independent   of  v  and  w. 

AAA,  /\ 

Proof.   Let  K  =   K     +   K„.   Then  K   is  contained  in  the  null  space  of  tt.  . 

Since  tt.  is  the  identity  mapping  on  L, 

(4.5.2)  K   n  L  =   3  . 

Now  consider  the  functional 

f(w,  v)  =  || w  -  v ||  „ 

L2(TR) 

defined  on  K   x  L.      Let  C  be  the  compact  subset 

{(w,  v):  ||w||       =  ||v||       =1} 
LZ(TR)      L^(TR) 

of  K   x  im      Since  f  is  continuous  it  attains  its  intimum  of  some  (w_,  v_)  G  c, 


but  by  (4.5.2)  we  have 


or 


|WQ  -  v  ||       >  c  >  0  , 
U  LZ(TR) 


|wQ  -  v  ||  2     >  c  1 1  w  1 1 

U    U  LZ(TR)       °  L^(TR) 


Since  f(w,  v)  is  minimal  at  (w  ,  v_) ,  (4.5.1a)  holds  on  all  of  C,  and 
then  on  K   x  L   by  linearity. 

For  (4.5.1b)  note  that  since  the  functions  in  K.   are  continuous 

1  ~  i  i 

and  piecewise  C  ,  the  only  functions  u  £  K.   for  which  |u|  ,  =0  are  the 

H 

constant  functions.   The  only  constant  function  in  K   is  the  function 

identically  zero,  since  K   is  contained  in  the  null  space  of  TT.  and  the 

function  values  at  the  vertices  of  T^  are  nodal  parameters  for  L.      Thus, 

K 
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• |  -  is  norm  on  K.      Because  of  this,  the  verification  of 
H 
(4.5.1b)  goes  through  exactly  the  same  as  for  (4.5.1a).  □ 


2 
For  any  triangle  T  C  IR  ,  each  affine  transformation 

F  :  T  ->■  T  induces  a  dual  mapping  F'  carrying  functions  on  T  back  to 
IK  T 

functions  on  T^ : 
R 

F^:  hV)  ->  H1^) 

by 

Fj(u)  =  u  o  FT  . 

One  can  then  ask  for  the  norm  of  this  dual  mapping  F'  either  as  a 

mapping  from  H  (T)  to  H  (T„)  or  from  L  (T)  to  L  (T0).   For  the  H 

seminorm  the  answer  is  that  |u|  ,     and  |u  °  F  |       are  essentially 

IT(T)  T  T1^) 

the  same  since  the  change  in  area  over  which  |  •  |  ..  is  computed  and 

H 
change  in  the  size  of  the  derivatives  are  opposite  effects  and  cancel 

perfectly  (in  two  dimensions).   Only  triangles  T  €E  T    are  considered 

in  the  following  lemma  both  because  they  are  the  only  ones  of  interest 

here  and  because  the  requirements  of  section  4.2  force  all  triangles 

in  T .  .  to  have  diameter  comparable  to  h.  , ,  simplifying  the  statement 
J-l  J-l 

of  the  lemma. 


Lemma  4.9.   Let  T  G  T   ,  .  j  >  2  and  let  F  be  an  affine  transformation 
J-l    - 

from  T_  onto  T.   If  u  e  H  (T) ,  then  u*  =  u  °  F  e  H  (TD)  and 
R  R 

(4.5.3a)         -  |u|      <  |u*|       <  c  |u|  . 

C    HX(T)       HX(TR)        HX(T) 
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(4.5.3b)  £  HI  2         iNll^Hz  1CHI2 

l/(T)     J      l/(TR)        LZ(T) 

for  fixed  c  >  1  independent  of  u,  j,  and  T. 

2 
Proof.   F  =  Ax  +  b  for  some  2x2  matrix  A  and  b  e  R  .   Then  as  shown 

in  Oden  and  Reddy  (1976,  eqs.  (6.167)  and  (6.168)),  for  any  s  >  0 
(A. 5.4a)       |u*|       <  ||a||S  Idet  A|~1/2  lul 

*  I       I    q  —  Hill  I  I     I    q  ■ 

HS(TR)  HS(T) 

(4.5.4b)      lul     <  Ha"1!!8  Idet  a|1/2  |u*| 

*  '  I   I   o       —  II       III  I        lie         ' 

HS(T)  HS(TR) 

where  ||A||  is  the  spectral  norm.   Now  let  dT  and  d   be  the  diameters  of 

R 

T  and  T  ,  respectively  and  let  p_  and  p_  be  the  diameters  of  the 
R  T      TR 

inscribed  circles  of  T  and  T  .   Then  as  given  in  Oden  and  Reddy, 
lemma  6.2,  the  following  simple  bounds  on  the  norms  of  A  and  A  "  apply 

llA||< 


PT 
R 


IKx|li 


We  also  have 


PT 


,       area  (T) 

det  A  = 


area  (TR) 


The  triangle  T  satisfies  the  angle  condition  (4.2.1)  uniformly  for 
T  S  t.  -  and  j  >_  2.   Using  this  condition  it  can  be  shown  by  simple 
geometric  arguments  that 
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—  d  <  p 

c   T  -  PT 

-  d   <  pm 
c   T  -  HT 

R     R 

^  d^  <  area(T)  <  c  d2 

\&\     <   area(TR)  <  c  d^   , 
R  R 

for  c  >  1  independent  of  T  £  T,    and  j  >_  2 .   Then  since  d   =  1  by 

R 
assumption  it  follows  that 

I|a||<  c  dT 

liA"1!!  <  e  d"1 

|det  A|1/2   <cdT 

|det  A|~1/2  fed"1. 
Combining  these  relations  with    (4.5.4)    yields    (4.5.3a)    and 

\\M,         <   dT  ||U*||  2  <e||n||  . 

l/(T)  L   (TR)  LZ(T) 

Since 

-  h.    ,    <   dm  <   c  h.    , 
c     j-1  -     T  -  j-1 

by  (4.2.  ),  and  p  =  2,  inequality  (4.5.3b)  follows  and  the  proof  is 

complete.  □ 

Lemma  4.9  above  shows  how  the  norm  of  a  function  u  on  a  triangle 

T  is  affected  by  transforming  it  to  a  function  u*  =  u  °  F^  on  T  . 

Lemma  4.8  showed  that  functions  in  K     or  K     were  not  well  approximated 

R     G 

by  those  in  L.      Combining  these  two  lemmas  we  should  be  able  to  show 
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ft     -1         -1 
that  for  any  triangle  T,  functions  in  K     °  F   or  K  °   F   are  not 

K    r      G    T 

well  approximated  by  those  in  L   °  F   .   This  is  equivalent  to  the 

statement  that  for  any  T  e  t    functions  in  M  |   cannot  be  well 

approximated  by  those  in  M.  .   ,  since  M,  |_  is  either  {0},  K     °  F   , 

j-1  1         j  1  R    T 

-1 
or  K     o  F   for  one  of  the  affine  transformations  F„,:  T_,  -*■   T.   Applying 
o     I  IK 

this  observation  to  all  triangles  T  e  T    yields  the  following  lemma, 
which  shows  that  functions  in  M._..  cannot  well  approximate  those  in  M.. 
This  lemma  was  shown  previously  by  Bank  and  Dupont. 


Lemma  4.10.   If  v  £  M.  -  and  w  £  M.  the  strengthened  Cauchy  inequaltiy 

(4.5.5)  a(v,   w)    <  y    |||v|||     |||w||| 

holds  for  fixed  y   £  (0,  1)  independent  of  v,  w  and  j  >_  2. 

Proof.   The  proof  will  be  done  locally  over  the  triangles  of  T .  _ . 

Let  T  S  t.  .  be  arbitrary.   Then  there  exists  an  affine  transformation 
J-1 

F  :  T  -*■  T.   Letting  w  €  $.  be  arbitrary  we  can  define  w*  w  °  F  .   Then 
T   R  2  *■ 

w*  £  K„  or  w*  £  K  even  when  T  £  T.  ,  and  w*  =  0.   Also  we  have 
R  6  j-1 

{v*  =  v  °  F  for  v  £  M.   }  c  L   where  the  inclusion  may  be  strict 

because  of  boundary  conditions  or  interelement  continuity  conditions. 

Thus, 

inf   |v  °  F  -  w*|       >  inf  |v*  -  w*|  1     >  c  |w*|  1 
v£M._1       T      HX(TR)   v*GL  HX(TR)  "       HX(TR) 

inf     ||v   o    F      -  w*||  ?  >     inf    ||v*  -  w*||  ?  >   C  ||w*||  , 

vSM._1  i  LZ(TR)        v*GL  L^(TR)  LZ(TR) 

where  the  inequalities  on  the  right  are  from  lemma  4.8.   By  lemma  4.9 
these  inequalities  can  be  converted  into 
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inf    |v  -  w|  1     ii  c  |w|  i 
v€M._         H  (T)         H  (T) 


inf 


v  -  w  II  „      >   c  ||w||  „     , 


^J-l        L2(T)   "      L2(T) 

where  the  constants  here  are  Independent  of  dT«   Combining  these 

relations  yields 

2  2  2 

inf   || v  -  w||  .     >_     inf    |v  -  w|  .     +  inf   ||v  -  w||  _ 

veM._1  h  (T)  ""  v£M.  h  (T)       vGM  L   (T) 

>   c(|w|2  +||w||2  )    =   c  ||w|| 

HX(T)  L    (T)  H    (T) 

Then  summing  over   the   triangles   of   x 

2  2 

inf       || v  -  w 1 1       =      inf  E  ||v  -  w||  , 

v€M.    ,  H  v£M.    ,    T^T.    ,  H    (T) 

3-1  j-1         j-1 

>  E  inf       ||v  -  w||2 

"  T^t.    ,    vSM.    .  H    (T) 

j-1  j-1 

>_  c  E       ||w||  =   c  ||w||  ,    . 

T^T  H^T)  H 

So  using  the  equivalence  of  the  H  and  energy  norms 

inf    |||v  -  w|||  >  c  MI  w  Ml  . 
J-1 
This  is  equivalent  to  (4.5.5).   To  see  this,  assume  |||w||  =  1.   Then 

inf        |||v  -  w|||    >      inf         |||v  -  w|||    >   c    |||w|||    =   c    , 

|||v|||=l  vSM 

vSM.    ,  J_i 

J-1 

so   if    |||v|||    =    |||w|||    =   1 

a(v,  w)    =   l-\   |||v-w|||2   <  1-   c-    (1-   c)    |||v|||     |||w|||    . 
The   general   case  of    (4.5.5)    follows  by   linearity.   □ 
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One  more  lemma  is  needed,  which  will  be,  for  our  purposes,  the 

most  important.   It  states  that  functions  in  M  ,  the  kernel  of  the 

interpolation  mapping  from  M   to  M,_1 ,  are  all  quite  oscillatory. 

Intuitively,  the  relatively  smooth  functions  in  M.  agree  quite  closely 

with  their  interpolate  in  M.  ,.   Only  the  most  oscillatory  functions  in 

M.  can  have  the  zero  function  as  their  interpolate.   This  observation 
J 

can  be  made  precise  quite  easily  through  the  use  of  lemma  4.9.   For  our 

purposes  a  function  is  "oscillatory"  when  its  H  norm  or  seminorm  is 

2  ~ 

much  larger  than  its  L  norm.   That  functions  in  M,  are  oscillatory  in 

this  sense  is  shown  by  the  following  lemma. 


Lemma  4.11.   If  w  €E  M .  ,  i  >  2  then 
j'   ~ 

(4.5.6)  -llwll  ,  <  h.    Iwl    -   < 

c  M     ii  2  -     3    i     II- 


c  ||w||  ?    , 
L"  J  L 


with  c  >  1  independent  of  w  and  j 


To 


Proof.   It  is  enough  to  show  this  locally  on  the  triangles  of  T .  - . 
see  this,  suppose  that  for  all  w  S  M.  and  T  £  T ._ 

(A. 5. 7)  i||w||  <  h      |w|  <   c  || w|| 

L  (T)     J     H  (T)        L  (T)  ' 

with  c  >  1  independent  of  w,  T,  and  j.   Then  for  any  w  £  M. 

\       Z    llwH22    1  hi    ^    I w|      <  c2   E    ||w||2     , 
c  tet      l  (t)    3   tgt._1   h  (T)     ^-i   L  (T) 

which  is  the  same  as  (4.5.6). 

To  show  (4.5.7),  let  w  G  M.  and  let  T  be  an  arbitrary  triangle 

in  T .  ..   Then  there  is  an  affine  transformation  F:  T_.  -*-   T  such  that 
J-l  R 

w*  =  w  o  F  is  in  K     or  K„.      Since   •   ,  is  a  norm  on  K_,  and  on  K     and 

R     G  tt1  K         b 

n 
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u  6  K     or  u  £   L 
R  6 


^  Hull  2  <l«l,  <e||n|| 

l/(TR)  H^CTjj)  l/(TR) 


with  c   >   1   independent  of  u.      Then 
1 


i|W*||  2  £    |w*|  <   c  ||w*||  2 

l/(TR)  HX(TR)  l/(TR) 


so  using  lemma  4.9 


-  II 


1  h,  lwl  i     1  c  llwll  2 
J     H  (T)         L 


L  (T)     J     H  (T)        L  (T) 
completing  the  proof.  □ 

4.6.   Approximat  ion 

In  this  section  the  approximation  inequality  (4.3.4)  needed 
to  justify  the  convergence  and  complexity  results  of  sections  4.3  and  4.4 
will  be  proven.   This  inequality  generalizes  a  result  of  Bank  and  Dupont 
for  quasi-uniform  grids.   Assume  the  elliptic  equation  is  H    regular 
for  a  ^  (0,  1).   Then  their  result  is  that  if  (M.}._   is  a  family  of 
quasi-uniform  spaces  and  S.  and  0.,    j  >_  2,  are  defined  as  before,  the 
following  inequality  holds  for  T)   €=  S .  , 


|||n-n|||£ce«/2  |||n|||   , 

where  r\   is  the  elliptic  projection  of  n.  on  M.  ...   A  domain  with  a  slit, 
the  worst  case  ordinarily  arising,  yields  a  =  -r-  and  their  result 
becomes 
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llln-n|||  <c  ej/4  |||.,|||  , 

the  same  as  the  result  here. 

The  result  here  generalizes  theirs  in  the  sense  that  it  applies 
to  locally  refined  grids  as  well  as  to  quasi-uniform  grids,  although  the 
exponent  of  6   is  usually  poorer  in  our  inequality.   This  is  of  secondary 
importance  for  multi-level  convergence  since  0(N)  complexity  can  be  shown 
for  any  positive  exponent  of  0..   Only  the  hidden  constant  in  the  0(N) 
bound  is  affected  by  the  exponent. 

Enough  has  been  said  about  this  approximation  result  that  no 
further  introduction  seems  necessary.   The  technique  of  proof  is  quite 
interesting  and  begins  with  a  sequence  of  strengthened  Cauchy  inequalities 
that  follow  more  or  less  directly  from  the  results  of  the  last  section. 
Each  of  these  inequalities  is,  in  effect,  a  lower  bound  on  the  dihedral 
angle  between  two  subspaces. 

Lemma  A .  12 .   Let  w  G  $ . ,  n.  G  S . ,  sGS.  ,  and  o  G  0 .    n  .   Then  the 
following  strengthened  Cauchy  inequalities  hold: 


(4.6.1a) 


a(n,  w)  <  c  QV2    |||n 


w    , 


(4.6.1b) 


a(s,  w)  <  c  Qlf_2±    |||s 


w 


(4.6.1c) 


a(n,  o)  :  c 


a/2  Q-i/2 


J      j-1 

Proof.   For  (4.6.1a)  let  r|  be  expanded  as 

JO.) 

3      (i) 

n=l 
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Then  the  corresponding   finite  element   data  is 


JO.) 
f  =       Z       a.   xfj)  ipfj)    , 


1-1       x     X 


with 


(4.6.2)  a(n,   40   =   (f,  <j>)    ,  <J>  e  Mj    . 


Then 


J(6.) 

i=l 


JO.) 


<  9.    A„J/        Z       at  X^J/  IU> 
.    .         ix        "    i 


J         1-1 
J(6.) 

=  e    \W     zJ    a2  |||^)|||2 
j      1-1 

-e    ^j)  IN"2  • 
j 


Thus 


(4.6.3)  ||f||<    (6     A^)1/2    IHnlH    . 

j 

Setting  (f>  =  w  in    (4.6.2) 

(4.6.4)  a(n,   w)   =    (f,   w)    <    ||f||||w||  <    (6.    ^j))1/2    |||n|||   ||w|| 

j 

so  using   lemma  4.11, 

(4.6.5)  a(n,   w)   <   c  h   (0      ^j))1/2    IHnlH     |||w||| 

<c  e)n  IHnlH   lllwlH  , 

using  the  eigenvalue  bound  (2.3.2).   The  proof  of  (4.6.1b)  is  the  same 
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(1-1)  Ni-1 

except  6    and  the  eigenfunctions  and  eigenvalues  {^.     }   .   and 

(1-1)  ^1-1 
{A  J    }  _   occur  throughout.   Thus,  instead  of  (4.6.5)  we  get 

a(s,  w)  <ch  (8    ^1))1'2  |||. HI  Hlwlll 

icp"1  i]H    lllslH    |||w||| 
ic6^   HI  s|||    i|j w|||    , 
since  by  assumption  p  =   2.      For    (4.6.1c)   we  have  as   in    (4.6.4) 

a(n,  o)  <  (e    xj  )1/2  IHnlll  ||o||  . 
j 

Expanding  o  in  eigenfunctions 

J_1        (1-1) 


Thus, 


N'-l 

|o|ll2=       ^        4  lll^lll2 
i=j(e._1)+i 

>e.  ,  tf-i)       V        ajll^ll2 

J  1-1 


=  e.  .  x.9"15  lloll2 


j-l     N.    , 

3-1 


so 
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a(n,  o)  <  (0    a^V/2  o       x^r1'2  Hlnlll   lllolH  . 
j  j-1 

By  the  eigenvalue  bound  (2.3.2) 

A<J>  <  c  hT2  . 
N.   -    J 

The  inequality  in  the  eigenvalue  bound  (2.3.2)  also  goes  the  other 

way.   One  way  to  see  this  is  to  expand  w  in  eigenf unctions  in  lemma  4.11: 

N. 
J 

w  =     I     a.   \b . 

i-1     X     X 
Then  the   left   side  of    (4.5.6)    is 

N.  N.  N. 

\\\l     a     !p.||<h     |   I     a    *   |        <  c  h     HI    Z     a     *   ||| 
i=l     x  J     i=l     1     x  „1  3       i=l 

rl 

where  the  last  inequality  is  from  the  equivalence  of  the  H  and  energy 


norms.   It  follows  easily  from  this  that 

A.""1*   >  e  h"2 


Yi  "     I'1  ' 


so  we   get 


a(n,  o)  <  c  p  ei/2  e"i(2  |||n|||   |||o||| 
<»ej/2e-^2  Hlnlll   li|o|||   , 

since  the  mesh  ratio,  p,  is  fixed.   This  completes  the  proof.  □ 


Each  of  equations  (4.6.1)  of  this  lemma  states  that  the  energy 
inner  product  of  a  relatively  smooth  function  with  a  relatively  oscillatory 
function  is  small  compared  to  the  product  of  their  norms.   In  other  words, 
the  angle  between  these  functions,  as  measured  by  the  energy  inner 


content 
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product,  is  large.   In  the  appropriate  limit  as  0  on  0.,  goes  to  zero 
each  of  equations  (4.6.1)  expresses  an  orthogonality  relation. 

Our  aim  is  to  show  that  the  function  in  S     can  be  well 

approximated  in  the  coarser  finite  element  space  M    when  0.  is  small. 

1 
Let  M,  ,  be  the  orthogonal  complement  of  M,  -  in  M,  with  respect  to  the 
j-l  j-l     J 

energy  inner  product.   We  approach  the  question  of  approximation  of 

functions  of  S.  in  M    by  the  back  door.   Instead  of  showing  directly 

that  M.  ,  contains  a  function  near  any  function  in  S . ,  we  show  that 

J-1  J 
1 

S.  and  M    are  nearly  orthogonal.  Since  S.  consists  of  relatively 

smooth  functions,  by  analogy  with  the  lemma  just  proven  we  should  try 

to  show  M.  ,  consists  of  oscillatory  functions.   This  is  the 

of  the  following  lemma. 


Lemma  4.13.   If  u  G  M,   ,  then 

u  =  s  +  o  +  w  , 

for  s  £  S .  , .  o  £  0 ,    ,,  and  w  £  M .  and 

J-l       J-l  J 

(4.6.7a)  |||s|||    <cQ][l    |||u|||    , 

(4.6.7b)  |||o|||    <   c    HI -u. HI    , 

(4.6.7c)  |||w|||    <   c    |||u|||    , 

with  c  >  0  independent  of  u  and  j. 

Proof.   Since  M.  -  M.  .,  (+)  M . ,  we  can  set 
J    J-l  v^   J 

U  =  V  +  w  , 

for  some  v  e  M.   ,  w  G  M..   Then  using  the  elementary  inequality 

xy  <  \   (x2  +  y2)  ,  x  ,  y  e  R  , 
and  lemma  4.5.3, 
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lllulll2-    l||v  +  W|||2=    |||v|||2+   |||w|||2  +  2a(v,  „) 

>  lllvlll2  +    ll|w|||2  -   2Y    |||v|||     |||w||| 

>  (1-Y  )(H|v|||2+  |||w|||2)    • 
Thus, 

(4.6.8a)  Hlvlll    <   c    |||u|||    , 

(4.6.8b)  '  |||v|||    <   c    |||u|||    , 

where  c  >  0  depends  only  on  the  constant  y  e  (0»  1)  from  lemma  4.10. 

1 
giving  (4.6.7c).   Now  since  u  =  v  +  w  £  M.  ..  , 

3-1 

a(v  +  w,  (())  =  0  ,  4>  G  M    » 
so  v  is  the  elliptic  projection  of  -w  on  M.  _.   Writing 


v  =  s  +  o  , 


for  s  £  S.  .  ,   o  €  0 .    ,  , 
J-l       J-l 


since 


a(s  +  w,  <$>)    =  0  ,  (J)  S  S.  -  , 
a(o  +  w,  <J>)  =  0  ,  (J)  G  0.    , 
S.    and  0.    are  energy  orthogonal.   That  is,  s  and  o  are  the 


elliptic  projections  of  -w  on  S.  .,  0 .      ,  respectively.   Then 

lllsHI2  =  -a(s,  w)  <  c  9^2  |||s|||  ilfwMI  , 
by  lemma  4.12,  (4.6.1b).   Thus,  using  (4.6.8b)  we  have 

lllslllf  c  6^    |||.|||   <c  eJ/J    |||u||| 

and 

III"  III  1  INII   <c    |||u|||    , 

completing   the  proof.  D 
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The  lemma  just  proved  shows  that  M._  ,  the  orthogonal 


complement  of  M.  n  in  M.,  is  contained  in 
J-l     J 


Vi  ©  Vi  ©  M: 


1 

and  for  any  u  £  M.  ..  ,  the  component  of  u  in  S.  ,  is  relatively  small 
J-l  J-l 

1 
when  6,  ,  is  small.   In  other  words,  M.  ,  consists  of  relatively 
J-l  J-l 

oscillatory  functions.   Since  S.  consists  of  smooth  functions  when 

J 

6.  is  small,  one  might  expect  that  the  inner  product  of  functions 

1 
in  S.  with  functions  in  M.  .  would  be  small  when  6.  was  small.  This 
J  J-l  J 

is  the  content  of  the  next  lemma. 


Lemma  4.14.   If  r\   £  S.  and  u  G  M.  n,  then 

J  J-l 

a(n,  u)  <  cej74  IHnlH   |||u|||  . 

Proof.   Since  0.  is  the  only  member  of  the  family  {6.}.  .,  of  parameters 

J  i  1=1 

occuring  either  implicitly  or  explicitly  in  the  statement  of  this 

lemma,  we  are  free  to  choose  the  remaining  parameters  {9.}.,.  arbitrarily. 

In  particular,  0.   E  (0,  1)  may  be  considered  a  free  parameter  until 

its  value  is  chosen  at  the  end  of  the  proof.   Let  r\   and  u  satisfy  the 

hypotheses.   According  to  lemma  4.13,  u  can  be  decomposed  as 

u  =  s  +  o  +  w  , 

with  s  e  S .  , ,  o  E  0 .  ,  and  w  e  M . .   Then 
J-l       J-l  J 

a(n,  u)  =  a(n,  s)  +  a(n,  o)  +  a(n,  w)  . 
Bounding  each  of  these  terms  separately 

a(n,  s)  <  IHnlll   |||s|||  <  c  eJ/J  |||n|||   |||u||| 

by  lemma  4.13,    (4.6.7a).      By  lemma  4.12,    (4.6.1c), 


Ill 


a(n,  o)  <  c  e1.'2  e-_x(2  IHnlll   |||o|||  <  e  ejfj  |||n|||   |||u|||   , 

where  the  last  inequality  is  from  lemma  4.13.   Similarly  by  lemma  4.12, 
(4. 6. Id), 

a(n,  w)  <  c  ej/2  llhlH   |||w|||  <c  q]/2  |||n|||   |||«|||  . 

Combining  these  inequalities 

a(n,  u)  <  c(e^  +  ei/2  9-1/2  +  ei/2,  |||nm   |||um   _ 

1/2 
Choosing  0.  n  =9.    completes  the  proof.  □ 
J-l    3 


Our  first  approximation  result  is  a  simple  consequence  of 
this  lemma. 

Theorem  4.15.   Let  r\  E.  S .    and  let  n  G  M.  ,  be  the  elliptic  proiection 
J  J-l 

of  ri  on  M.  ..  , 
3-1 

a(n  -  n,  *)  =  0  ,  cf)  e  M    . 
Then 

llln-nll!  £cej/4  |||n||,   . 

Proof.      We  have 

|||n-  n|||2  -  a(n,  n  -  n)  <  c  e^4  |||n|||    |||n  -  n|||   , 

by  lemma  4.14,  since  r\   -  n  G  M .   .   Dividing  by   |  ri  -  n  |   the 
result  follows.  □ 

An  important  feature  of  the  result  just  proven  is  that  the 
constant  c  occurring  here  depends  entirely  on  local  properties  of  the 
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finite  element  spaces  M    and  M  .   In  fact,  c  depends  only  on  the 

constants  in  the  eigenvalue  bounds  and  on  those  of  the  last  section 

on  interpolation.   As  a  consequence,  it  is  readily  computable. 

Moreover,  the  approximation  property  just  shown  does  not  depend 

significantly  on  either  the  domain  or  the  boundary  conditions.   This 

is  in  contrast  to  the  result  of  Bank  and  Dupont  where  both  the 

exponent  of  8.  and  the  constant  occurring  depend  on  the  regularity 

of  the  PDE  and  hence  on  the  domain. 

The  approximation  result  just  proven  is  almost  the  result 

we  needed  in  section  4.3.   Its  defect  is  that  it  is  expressed  in 

N. 
terms  of  the  subspace  5.  and  indirectly  the  eigenfunctions  ity  .}      ., 

J.  5  x 

rather  than  in  terms  of  S.  and  indirectly  0i>.}.  , .   Since  it  would 

J  i  i=l 

be  difficult  to  analyze  the  simultaneous  displacement  iteration  of 

N. 
section  4.3  in  terms  of  the  eigenfunctions  {i|0._,  we  must  prove 

the  analog  (4.3.4)  of  theorem  4.15  which  has  S.  replaced  by  S..   Exactly 

the  same  thing  occurred  in  chapter  three.   While  it  may  be  possible 

to  get  the  required  analog  directly  from  theorem  4.15,  it  is  easier  to 

repeat  the  steps  in  its  proof  with  slight  modifications.   Actually, 

this  could  have  been  done  in  the  first  place,  but  it  seemed  more 

natural  this  way.   First  we  show  the  analog  of  lemma  4.12.   Only  S. 

need  be  replaced  by  S . .   There  is  no  necessity  to  replace  S.  -  and 

0.    1  since  these  subspaces  play  only  an  intermediate  role. 

Lemma  4.16.   Let  wGM.,r)eS.,.sGS.  .  and  o  G  0 .  ..   Then  the 
J       J       3-1         J"1 

following  strengthened  Cauchy  inequalities  hold: 
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(4.6.9a)  a(n,   w)    <   c  ej/2    |||n|||     |||w|||    , 

(4.6.9b)  a(s,   w)   <   c   8^2    |||s|||     |||w|||     , 

(4.6.9c)  aOn,    o)   <   c  Q1.'2  6^{2    |||n|||     |||o|||    . 

Proof.   Equation  (4.6.8b)  is  the  same  as  (4.6.1b)  so  is  already  established, 
For  (4.6.8a)  expand  n  in  the  eigenfunctions  {$.   }, 

JO.) 

J  ~(i) 

n  =     E         a.  <KJ'    , 

i=l         X     X 
and  let   f  be  given  by 

JO.) 

f  =       I       a.   X?j)   ifj)    . 
i-1       x     x         x 

Then 

(4.6.10)  a(n,    <$>)   -  b(f,    (j>)    ,    <J>  e  M.    • 

Now  instead  of  (4.6.3)  we  get  by  a  similar  argument 

IUIU  ce1  ^j))1/2  HlnlH   , 
j 

so  with  f  =  w  in    (4.6.10), 

(4.6.11)  a(n,   w)   =  b(f,   w)    <    (9.    X^j))1/2    |||n|||   ||w|^    . 

Then  as  in  (3.3.3),  since  the  eigenvalues  of  the  mass  matrix  lie  on 
an  interval  [ — ,  o]    for  cr  >  1  independent  of  j  , 

(4.6.12)  ||»H,  <  O1'2  ||„||  <  e  01/2  hj    IH-lll     , 

where  we  used  lemma  4.11.   Combining  (4.6.11)  and  (4.6.12)  and  using 
the  fact  that  a  is  fixed, 
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a(n,  w)  <ch  (9    ^>)1/2  Ulnlll   Hl.lll 

j 

<c  e]n  lllnllllllwlll 

by  the  eigenvalue  bound  (2.3.2).   For  (A. 6. 9c)  we  have  as  in  (A. 6. 11) 

a(n,   o)   <    (6     tfh1/2    HlnlH   ||o|^  <  a1/2    (9     A^)1/2    |||n|||   ||o  ||  . 
j  J        j 

Then  as   in    (4.6.6) 

NI<<Vi^d~1>ri/2  INolll  , 

3       j-i 


so 


a(n,o)<a1/2(e1x«')1/2  (6..JO)  ri 

j  j-l 


HlnlH    lllo 


<c(6.  h"2)172  (9.^  h-2_±)-1/2  HlnlH    HlolH   , 

using  the  eigenvalue  bounds  as  before.   Then  since  the  mesh  ratio   is 
fixed, 

a(n,  .)ic8J,!  e-i(2  HlnlH   |||o||| 

as  required.  □ 


Since  we  are  retaining  S.  -  and  0     .    and  only  replacing  5. 

by  S.,  there  is  no  need  to  alter  lemma  4.13.   Then  following  the  same 
J 

proof  as  for  lemma  4.14,  but  using  the  strengthened  Cauchy  inequalities 
of  lemma  4.16  rather  than  lemma  4.12,  we  get  the  analog  of  lemma  4.14. 
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Lemma  A.  17.      If  ri  e  S  and  u  £  M.    ,    then 
J-l 

a(n,  u)  <  c  Qx./U   HlnlH    |||u|||    . 

The  required  approximation  theorem  follows  easily  from  this 
lemma  just  as  the  analogous  theorem  4.15  followed  from  lemma  4.14. 

Theorem  4.18.   Let  J]   £  S.  and  let  rj  £  M.  ..  be  the  elliptic  projection 

of   T]   on  M.       , 

a(n  -  n,  <}))  =  o  ,  <f>  e  m._1  . 

Then 

(4.6.13)  Illn-n|||  <«  e^4  |||n|||  . 

With  this  approximation  result,  the  convergence  and  complexity 
results  of  sections  4.3  and  4.4  are  now  completely  justified.   As  for 
theorem  4.15,  the  constant  c  occurring  here  depends  only  on  local 
properties  of  the  finite  element  spaces.   For  particular  spaces  one 
could  easily  compute  this  constant  and  determine  a  rigorous  upper  bound 
on  the  rate  of  convergence  of  the  multi-level  iterations  of  section  4.3. 
This  generalizes  the  heuristic  local  Fourier  mode  analysis  to  irregular 
finite  element  grids,  and  is,  moreover,  completely  rigorous.   The  hidden 
constants  in  the  0(N)  complexity  bounds  of  section  4.4  are  not  so 
easily  computed  since  they  depend  on  the  generally  unknown  constant  in 
the  error  estimate  (4.4.1).   However,  for  most  purposes  it  is  enough 
to  have  a  bound  on  the  spectral  radius  of  the  multi-level  iteration,  which 
this  analysis  gives  since  this  is  what  determines  an  iteration's 
practicality. 
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